B  and  F  Projection  Methods  for  Nearly 
Incompressible  Linear  and  Nonlinear 
Elasticity  and  Plasticity  using  Higher-order 
NURBS  Elements 

T.  Elguedj 1  ,Y.  Bazilevs1  ,V.M.  Calo1  ,T.J.R.  Hughes2 

Institute  for  Computational  Engineering  and  Sciences,  The  University  of  Texas  at 
Austin,  201  East  2fth  street,  1  University  Station  C0200,  Austin,  TX  78712-0027, 

USA 


Abstract 

This  paper  presents  projection  methods  to  treat  the  incompressibility  constraint 
in  small  and  large  deformation  elasticity  and  plasticity  within  the  framework  of 
Isogeometric  Analysis.  After  reviewing  some  fundamentals  of  isogeometric  analysis, 
we  investigate  the  use  of  higher-order  Non-Uniform  Rational  B-Splines  (NURBS) 
within  the  B  projection  method.  The  higher-continuity  property  of  such  functions 
is  explored  in  nearly  incompressible  applications  and  shown  to  produce  accurate 
and  robust  results.  A  new  nonlinear  F  projection  method,  based  on  a  modified 
minimum  potential  energy  principle  and  inspired  by  the  B  method  is  proposed 
for  the  large-deformation  case.  It  leads  to  a  symmetric  formulation  for  which  the 
consistent  linearized  operator  for  fully  nonlinear  elasticity  is  derived  and  used  in  a 
Newton-Raphson  iterative  procedure.  The  performance  of  the  methods  is  assessed 
on  several  numerical  examples,  and  results  obtained  are  shown  to  compare  favorably 
with  other  published  techniques. 
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1  Introduction 


1.1  Preamble 


Nonlinear  finite  element  structural  analysis  is  dominated  by  the  use  of  low- 
order  “displacement”  elements  that  are  specially  designed  to  avoid  volumet¬ 
ric  or  incompressible  locking.  Deformations  that  involve  very  small  volume 
changes  occur  for  rubber-like  materials,  and  the  clastic-plastic  response  of 
metals  and  undrained  soils.  Elastic-plastic  analysis  is  ubiquitous  in  many  en¬ 
gineering  disciplines  and,  thus,  it  is  essential  that  finite  elements  are  able  to 
accurately  represent  nearly  incompressible  deformations.  Standard  elements 
have  difficulty  or  fail  entirely  in  this  situation.  A  full  discussion  of  the  early 
experiences  and  development  of  effective  methodologies  is  presented  in  Hughes 
[29],  Chapter  4.  It  would  seem  at  this  stage  of  the  development  of  finite  cle¬ 
ment  technology  that  higher-order  approaches  would  play  an  important  role 
in  nonlinear  structural  mechanics,  but  this  is  not  the  case.  The  only  approach 
that  claims  any  success  is  the  p- method,  in  which  the  polynomial  order  within 
elements  is  increased  on  a  fixed  mesh  (see  [62,  63]).  The  claim  for  standard 
higher-order  C°-continuous  finite  elements  is  that  volumetric  locking  is  allevi¬ 
ated  as  the  element  polynomial  order  is  increased.  This  seems  to  be  the  case, 
but  there  is  evidence  that  the  accuracy  of  the  solution  at  any  fixed  polyno¬ 
mial  order  is  far  from  optimal  (some  confirmation  of  this  observation  will  be 
presented  herein).  In  addition,  numerical  experience  indicates  standard  higher- 
order  elements  are  much  more  “fragile”  than  low-order  elements.  This  lack  of 
robustness  is  particularly  apparent  in  nonlinear  dynamic  analysis  of  structures 
involving  contact  and  impact,  and  subject  to  high  wave  number  inputs,  such 
as  blast  waves.  The  reasons  for  this  have  never  been  fully  explained.  How¬ 
ever,  recent  vibration  and  wave  propagation  studies  (see  Cottrell  et  al.  [18] 
and  Hughes  et  al.  [37])  have  revealed  that  the  higher  modes  produced  by  the 
p-method  diverge  with  p.  That  is,  whereas  formal  accuracy  is  increased,  the 
improvement  is  confined  to  lower  modes,  while  at  the  same  time  the  higher 
modes  get  worse  as  p  increases.  This  may  explain  why  robustness  decreases 
with  p. 

An  alternative  higher-order  approach  has  recently  emerged  from  Isogeometric 
Analysis  (see  Hughes  et  al.  [32]),  namely  /;- refinement,  in  which  discretiza¬ 
tions  of  order  p  achieve  Cp_1  continuity  on  “patches”  (roughly  speaking,  sub- 
domains).  It  has  been  shown  that  fc-rehned  meshes  behave  entirely  differently 
than  standard  finite  element  methods  with  respect  to  higher  modal  compo¬ 
nents.  In  fact,  in  some  cases  all  discrete  modes  converge  to  exact  ones  and 
nearly  spectral  accuracy  is  achieved.  To  us,  this  suggests  that  robust  and 
higher-order  accurate  finite  element  methods  applicable  to  nonlinear  struc¬ 
tural  analysis  may  be  a  possibility,  and  this  study  represents  a  first-step  in 
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this  direction. 


At  the  same  time,  we  need  to  deal  with  the  locking  problem  in  both  small-  and 
finite-deformation  regimes.  We  feel  the  most  widely  applicable  and  practically 
useful  approach  is  by  way  of  a  pure  displacement  formulation  (he.,  no  pressure 
degrees  of  freedom)  as  is  generally  employed  in  large-scale  structural  analysis 
programs  (see  for  example  [43,  44]).  In  order  to  achieve  good  behavior,  we 
feel  it  is  imperative  to  use  some  form  of  “projection1'  to  reduce  the  number  of 
volumetric  constraints.  This  is  absolutely  essential  for  lower-order  elements, 
and  very  important  for  higher-order  elements  as  well  (as  we  will  show  herein). 
Designing  projection  schemes  for  higher-order  elements,  as  far  as  we  are  aware, 
has  not  been  pursued  previously  in  the  literature.  The  B  scheme  (see  Hughes 
[28])  is  a  formalism  that  utilizes  projection  and  we  propose  a  new  family  of 
higher-order  B  schemes  in  the  sequel.  Engineering  and  mathematical  points  of 
view  are  both  relevant  to  this  pursuit,  but  are  not  entirely  in  accord.  In  these 
cases,  we  favor  the  view  gained  from  engineering  experience  (de  gestibus  non 
est  disputandum) .  In  particular,  some  elements  produced  by  our  scheme  are 
definitely  known  not  to  satisfy  the  well-known  Ladyzenskaya-Babuska-Brezzi 
(LBB)  stability  condition  (see,  e.g.,  [10,  11,  62]),  but  nevertheless  are  favored 
and  widely  utilized  in  engineering  applications,  the  prime  example  being  the 
mean  dilatation,  bilinear  quadrilateral  element  “Qi/Qo”  ■  It  should  be  men¬ 
tioned  here  that  this  element  fails  to  satisfy  the  LBB  condition  only  “slightly” 
so  is  not  altogether  mathematically  abhorrent.  We  also  wish  to  emphasize  we 
are  certainly  not  opposed  to  satisfying  the  LBB  condition,  and  some  of  our 
other  constructs  may,  at  least  under  certain  situations,  but  we  also  want  to 
develop  methods  that  have  the  highest  probability  of  being  used  in  engineering 
applications.  This  is  the  key  issue  in  the  approaches  we  have  proposed.  The 
large-deformation  counterpart  involves  projection  of  the  deformation  gradi¬ 
ent,  a  so-called  F-scheme,  involving  a  product  decomposition  into  volumetric 
and  deviatoric  factors.  There  has  been  only  limited  use  of  such  schemes  in 
the  literature,  and  in  our  opinion  none  is  completely  satisfactory.  More  about 
this  later.  A  second  contribution  of  this  work  is  the  development  of  a  new 
F  formulation  based  on  a  modified  minimum  potential  energy  principle.  We 
derive  the  variational  equations  and  a  consistent  tangent  operator  for  gen¬ 
eral  hyperelastic  materials.  There  are  a  number  of  “unexpected”  terms  in  the 
consistent  tangent  and  it  is  symmetric,  a  theoretical  improvement  over  some 
previous  F  formulations.  We  also  perform  a  number  of  numerical  calculations 
on  quasi-static,  small-  and  large-deformation  test  problems.  The  calculations 
are  very  encouraging  and  seem  to  validate  the  new  methods. 
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1.2  Literature  survey 


Several  techniques  are  available  in  the  literature  to  deal  with  incompressibil¬ 
ity.  The  mixed  Galerkin  finite  element  formulation  is  one  of  the  most  com¬ 
mon.  However,  the  order  of  approximation  for  both  pressure  and  displace¬ 
ment /velocity  variables  cannot  be  chosen  arbitrarily.  The  method  must  ful¬ 
fill  the  LBB  condition  to  ensure  stability  and  optimal  convergence  (see,  e.g., 
[10,  11,  62])  or  must  be  used  in  conjunction  with  stabilization  techniques  which 
allow  the  use  of  a  wide  variety  of  interpolations,  including  equal-order.  These 
stabilization  techniques  have  been  utilized  with  success  in  the  past  twenty 
years  in  fluid  mechanics;  the  most  commonly  used  are  Streamline  Upwind 
Petrov-Galerkin  (SUPG,  see,  e.g.,  [12,  33])  and  Galerkin  Least  Squares  (GLS, 
see,  e.g.,  [34]).  In  more  recent  work,  these  ideas  have  been  successfully  ex¬ 
tended  to  nonlinear  solid  mechanics  [42,  46,  49].  Another  possible  form  of 
stabilization  for  displacement  based  formulations  is  “hourglass”  control  used 
in  conjunction  with  reduced  integration  rules  (see,  for  example,  [9,  25,  50,  51]). 

A  different  approach  to  handle  the  incompressibility  condition  is  the  incom¬ 
patible  mode  technique  proposed  by  Wilson  et  al.  [65],  and  later  generalized 
as  the  enhanced  or  assumed  strain  method  by  Sirno  et  al.  [57,  59]  within  a 
three-field  variational  principle.  Although  this  method  has  exhibited  hourglass 
instabilities  in  some  cases,  several  improvements  have  been  proposed  to  over¬ 
come  these  effects,  in  particular,  by  using  an  alternative  enhancement  strategy 
[2],  or  a  mixed-enhanced  strain  formulation  [40,  41], 

The  approach  of  greatest  engineering  interest  is  to  use  a  simple  pure- 
displacement  formulation.  This  idea  goes  back  in  the  1970’s  with  the 
so-called  “reduced”  and  “selective”  integration  methods,  and  B  technique 
[29],  for  which  equivalence  theorems  with  mixed  methods  have  been  obtained 
[28,  35,  45].  If  seen  in  a  more  general  way,  as  a  strain  projection  technique 
[29,  30],  one  can  include  within  the  same  B  framework  the  well  known 
mean-dilatational  formulation  of  Nagtegaal  et  al.  [47].  This  formulation  is 
widely  used  in  large-scale  and  commercial  codes  [43,  44],  The  major  numerical 
developments  on  F-type  methods,  the  generalization  of  B-type  methods  to 
hyperelastic  finite-deformation  formulation,  have  been  done  by  Hughes  et  al. 
[38],  Simo  et  al.  [60]  (within  a  three-field  Hu-Washizu  variational  principle), 
and  de  Souza  Neto  et  al.  [20,  22],  who  proposed  an  F  technique  for  linear 
elements.  These  works  exhibit  some  shortcomings:  a  consistent  tangent  was 
not  derived  in  [38],  it  was  never  published  in  the  open  litterature,  and  not 
pursued  further;  [60]  is  a  mixed  method  requiring  pressure  unknowns;  and 
[20,  22]  leads  to  a  non-symmetric  consistent  tangent. 

The  development  of  isogeometric  analysis,  recently  proposed  by  Hughes  et 
al.  [32],  introduced  a  new  numerical  method  in  an  attempt  to  improve  ge- 
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ometry,  solution  representation  and  mesh  refinement  compared  with  finite 
element  analysis.  Isogeometric  analysis  is  based  on  the  geometric  primitives  of 
computer  graphics  and  Computer-Aided  Design  (CAD).  The  first  implemen¬ 
tation  of  isogeometric  analysis  was  based  on  Non-Uniform  Rational  B-Splines 
(NURBS),  a  standard  technology  in  CAD  [15,  24,  48,  52],  It  was  first  applied 
to  three  dimensional  linear  elasticity  and  advection-diffusion  problems  in  [32] . 
The  ability  of  NURBS  to  represent  precise  geometry  was  successfully  applied 
to  solve  problems  of  vascular  fluid-structure  interaction  with  patient-specific 
geometries  in  [6,  66].  The  mathematical  theory  of  NURBS  based  isogeometric 
analysis  was  begun  in  [7],  where  stability  proofs  and  optimal  error  estimates 
were  obtained  assuming  /j-refinement.  An  interesting  feature  of  isogeometric 
analysis,  not  shared  by  finite  element  analysis,  is  so-called  /^-refinement,  which 
simultaneously  increases  polynomial  order  and  continuity  of  the  basis  func¬ 
tions.  Using  this  feature,  isogeometric  analysis  attains  better  accuracy  than 
finite  elements  for  structural  vibrations  [17,  18].  The  higher  continuity,  be¬ 
yond  the  geometric  interest,  has  also  proved  advantageous  for  turbulent  flows 
problems  [1].  The  above  developments  and  observations  suggest  NURBS-based 
isogeometric  analysis  may  offer  some  new  possibilities  in  large  scale  structural 
mechanics  applications. 


1.3  Outline 


The  paper  is  organized  as  follows:  In  Section  2  the  basics  of  isogeometric  anal¬ 
ysis  are  recalled  with  an  emphasis  on  k- refinement.  In  Section  3,  the  use  of 
high-order  NURBS  within  a  projection  technique  is  studied  in  the  geometri¬ 
cally  linear  case  with  a  B  method  to  investigate  the  choice  of  approximation 
and  projection  spaces  with  NURBS.  The  performance  of  the  proposed  strat¬ 
egy  is  demonstrated  on  well-known  numerical  examples.  Section  4  presents  the 
new  F  formulation.  The  consistent  tangent  operator  of  the  Newton- Raphson 
iterative  scheme  for  the  fully  nonlinear  elastic  case  is  derived  in  detail  in  Ap¬ 
pendix  B.  The  performance  of  the  new  method  is  compared  with  some  of  the 
techniques  cited  earlier  on  several  numerical  examples  in  Section  5,  and  the 
good  behavior  of  the  new  methodology  is  again  observed.  Finally,  Section  6 
presents  some  conclusions  and  possible  directions  of  future  work. 


2  Brief  overview  of  Isogeometric  Analysis 

This  section  gives  a  brief  overview  of  NURBS-based  isogeometric  analysis.  A 
more  detailed  description  can  be  found  in  [17,  32],  NURBS  functions  are  a 
standard  tool  in  CAD  to  describe  and  model  curves  and  surfaces  [15,  24,  48, 
52], 
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2.1  B-Splines 


NURBS  are  built  from  B-Splines.  B-Splines  are  piecewise  polynomial  curves 
composed  of  a  linear  combination  of  B-Spline  basis  functions.  The  B-Spline 
parametric  space  is  local  to  “patches”  in  contrast  with  finite  elements  in  which 
each  element  carries  its  own  parametrization.  Patches  may  be  thought  of  as 
subdomains. 


2.1.1  Knot  Vectors 

A  knot  vector  in  one  dimension  is  a  set  of  coordinates  in  the  parametric 
space  written 

1=1  =  {‘Cl)  •••)  Cn+p+l})  (1) 

where  ^  6  1  is  the  ith  knot,  i  is  the  knot  index,  i  =  l,2,...,n  +  p+l,pis  the 
polynomial  order  of  the  B-Spline  and  n  is  the  number  of  basis  functions  cor¬ 
responding  to  it.  The  knots  partition  the  parameter  space  into  elements,  and 
the  interval  [£i,  £n+p+i]  constitutes  a  patch.  A  knot  vector  is  said  to  be  uni¬ 
form  if  its  knots  are  uniformly  spaced  and  non-uniform  otherwise.  Knot 
values  may  be  repeated,  that  is,  more  than  one  knot  may  take  a  given  value. 
The  multiplicities  of  knots  have  important  implications  on  the  continuity  of 
the  associated  B-Spline  functions.  A  knot  vector  is  said  to  be  open  if  its  first 
and  last  knots  are  repeated  p  +  1  times.  In  what  follows,  we  always  employ 
open  knot  vectors.  Basis  functions  formed  from  open  knot  vectors  are  inter- 
polatory  at  the  ends  of  the  parametric  interval  [£l,  £n+p+i],  but  are  not,  in 
general,  interpolatory  at  the  interior  knots.  Open  knot  vectors  enable  patches 
to  be  assembled  in  essentially  the  same  way  elements  are  assembled  in  finite 
element  analysis. 


2.1.2  Basis  functions 


B-Spline  basis  functions  are  defined  recursively  starting  with  piecewise  con¬ 
stants  (p  —  0);  given  a  knot  vector  S 


Ni,  O(0  = 


if  &<£<&+i> 

0  otherwise. 


For  p  >  1,  they  are  defined  by: 


V,PK)  =  -AJyiv^g)  +  W1  (  iv-i+1J_l(o. 


£i+p  & 
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i+'P+ 1 


-6 


i+ 1 


(2) 

(3) 


An  example  for  a  uniform  knot  vector  is  presented  in  Figure  1.  Note  that  for 
p  —  0  and  p  —  1,  the  basis  functions  are  the  same  as  for  standard  piecewise 
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constant  and  linear  finite  element  functions. 


Fig.  1.  Basis  functions  of  order  0,1,2  for  an  open,  uniform  knot  vector 
S  =  {0,1, 2,  3, 4,...}. 


An  example  of  quadratic  basis  function  for  an  open  non-uniform  knot  vector 
is  presented  in  Figure  2.  Note  that  the  basis  functions  are  interpolatory  at 
the  ends  of  the  interval  and  also  at  £  =  4,  the  location  of  a  repeated  knot, 
where  only  C°-continuity  is  attained.  Basis  functions  of  order  p  have  p  —  m* 
continuous  derivatives  across  knot  £$,  where  m;  is  the  multiplicity  of  the  value 
of  in  the  knot  vector.  When  the  multiplicity  of  a  knot  is  exactly  p,  the 
basis  is  interpolatory  at  that  knot.  When  the  multiplicity  is  p+  1,  the  basis  is 
discontinuous  and  the  patch  is  split  into  two  separate  patches. 


Fig.  2.  Quadratic  basis  functions  for  an  open,  non-uniform  knot  vector 
£  =  {0,0, 0,1, 2, 3, 4, 4, 5, 5, 5} 

An  important  property  of  B-Spline  basis  functions  formed  from  an  open  knot 
vector  is  that  they  constitute  a  partition  of  unity,  that  is,  V£: 

n 

EAri,r«)  =  l-  (4) 

i= 1 

This  feature  is  also  shared  with  finite  elements  and  meshless  methods.  Another 
interesting  property  is  that  the  support  of  each  Ni  p  is  compact  and  contained 
in  the  interval  [£i,£i+p+i].  Finally,  one  can  note  that  each  basis  function  is 
point-wise  non-negative  over  the  entire  domain  of  definition:  >  0  V£.  This 

means  that  all  of  the  entries  of  the  mass  matrix  will  be  positive,  which  has 
implications  for  developing  mass  lumping  schemes  (see  [18]  for  an  initiatory 
investigation). 
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2.1.3  B- Spline  curves 


B-Spline  curves  in  are  constructed  by  taking  linear  combinations  of  B- 
Spline  basis  functions.  The  vector-valued  coefficients  of  the  basis  functions 
are  referred  to  as  control  points.  They  are  analogous  to  nodal  coordinates 
in  finite  element  analysis  in  that  they  are  the  coefficients  of  the  basis  functions. 
However,  the  non-interpolatory  nature  of  the  basis  does  not  lead  to  the  usual 
finite  element  geometric  interpretation  of  the  control  points.  The  piecewise 
linear  interpolation  of  the  control  points  defines  the  control  net.  Given  n 
basis  functions,  Ni>p,  i  =  1,2 ,  ...,n,  and  n  corresponding  control  points  G 
Md,  i  —  1,2,  ...,n  a  piecewise  polynomial  B-Spline  curve  is  given  by: 

n 

C(  0  =  E  mAOBi.  (5) 

i= 1 

An  example  is  shown  in  Figure  3  for  the  quadratic  basis  illustrated  in  Fig¬ 
ure  2.  Note  that  the  curve  is  interpolatory  at  the  first  and  last  control  points, 
a  general  feature  of  a  curve  built  from  an  open  knot  vector.  The  curve  is 
also  interpolatory  at  the  sixth  control  point.  This  is  due  to  the  fact  that  the 
multiplicity  of  the  knot  £  =  4  is  equal  to  the  polynomial  order  there,  that  is, 
m,5  =  p  =  2.  Note  also  that  the  curve  is  tangent  to  the  control  polygone  at  the 
first,  sixth  and  last  control  points.  The  curve  is  Cp~l  =  C1-continuous  every¬ 
where  except  at  the  repeated  knot  £  =  4  where  it  is  Cp~ms  =  C°-continuous. 


(a)  Curve  and  control  points  (b)  Curve  and  mesh  denoted  by  knot  locations 

Fig.  3.  Piecewise  quadratic  B-Spline  curve  in  Mrf.  a)  Control  point  locations  are 
denoted  by  #’s.  b)  Knots,  which  define  the  mesh  by  partitioning  the  curve  into 
elements,  are  denoted  by  B’s. 


The  properties  of  B-Splinc  curves  follow  directly  from  the  properties  of  their 
basis  functions.  For  example,  B-Spline  curves  have  continuous  derivatives 
through  order  p  —  1  in  the  absence  of  repeated  knots  or  control  points.  Re¬ 
peating  a  control  point  or  a  knot  k  times  decreases  the  number  of  continuous 
derivatives  by  k.  An  affine  transformation  of  a  B-Spline  curve  is  obtained 
by  applying  the  transformation  to  the  control  points.  This  turns  out  to  be 


the  essential  property  to  guarantee  satisfaction  of  so-called  “patch-tests”,  as 
discussed  in  [32],  This  property  is  referred  to  as  affine  covariance. 


2.1.4  h-  and  p-refinement:  knot  insertion  and  order  elevation. 

The  analogue  of  ^-refinement  is  knot  insertion.  Knots  can  be  inserted  with¬ 
out  changing  the  curve  parametrically  or  geometrically.  Given  a  knot  vector 
H  =  {£1, £n+p+i},  let  £  G  [£fc,£fc+i[  be  a  desired  new  knot.  The  new  n  +  1 
basis  functions  are  formed  recursively,  using  Eqs.  (2)  and  (3),  with  the  new 
knot  vector  S  =  {£1,  •••,  £fc,  £,  £fc+i,  •••,  £n+p+i}.  The  new  n  +  1  control  points 
{Bi,  ...,Bn+ 1}  are  formed  from  the  original  control  points  {£>1,  by 


Bi  ociBi  -(-  (1  1, 


(6) 


where 


1 


0 


if  1  <  i  <  k  —  p, 

if  k  —  p  +  1  <  i  <  k, 
if  k+l<i<n  +  p  +  2. 


(7) 


Knot  values  already  present  in  the  knot  vector  can  be  repeated.  However,  this 
reduces  the  continuity  of  the  basis  at  the  corresponding  knot.  The  continuity 
of  the  curve  is  preserved  by  choosing  the  control  points  using  Eqs.  (6)  and  (7). 


/i-refinement  p-refinement 


knot 

vector 


-={0,0,0, 1,1,1} 


H  =  {0,0, 0,0.5, 1,1,1} 


H  =  {0, 0,0, 0,1, 1,1,1} 


basis 

functions 


Fig.  4.  Curve,  control  points,  knot  vectors  and  basis  functions  after  knot  insertion 
and  order  elevation  of  a  quadratic  curve.  Control  points  are  denoted  by  #’s. 
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The  mechanism  for  implementing  /^-refinement  is  order  elevation.  The  poly¬ 
nomial  order  of  basis  functions  can  be  increased  without  changing  the  curve 
parametrically  or  geometrically.  Note  that  each  unique  knot  value  in  H  must 
be  repeated  in  order  to  preserve  discontinuities  in  the  pth  derivatives  of  the 
curve  that  is  being  elevated.  The  number  of  new  control  points  depends  on 
the  multiplicities  of  existing  knots.  As  for  knot  insertion,  the  solution  space 
spanned  by  the  elevated  basis  contains  the  space  spanned  by  the  original  ba¬ 
sis.  Thus,  it  is  possible  to  order  elevate  without  changing  the  geometry  and 
the  parametrization  of  the  B-Splinc  curve  (see  [48,  52]  for  further  details  and 
algorithms) . 

An  example  of  h-  and  p- refinement  is  presented  on  Figure  4.  The  original 
curve  consists  of  quadratic  B-Splines  with  knot  vector  S  =  {0,  0,  0, 1, 1, 1}. 
The  original  curve,  knot  vector  and  basis  functions  are  shown  on  the  left. 
A  new  knot  is  inserted  at  £  =  0.5.  The  new  curve,  knot  vector  and  basis 
functions  are  shown  in  the  center  of  the  figure.  Note  that  the  curve  remains 
unchanged  but  that  basis  functions  and  control  points  are  different;  there  is 
one  more  of  each.  Next,  the  curve  is  order  elevated  once;  the  new  curve,  control 
points  and  basis  functions  can  be  seen  on  the  right  of  the  figure.  This  time 
the  multiplicity  of  the  knots  is  increased  by  one.  The  location  and  number  of 
control  points  change,  but  the  curve  remains  the  same.  There  are  now  four 
cubic  basis  functions.  Despite  the  same  number  of  basis  functions,  note  that 
the  locations  of  the  new  control  points  obtained  from  knot  insertion  and  order 
elevation  are  different. 


2.1.5  k-reftnement:  higher- order  and  higher  continuity 

The  last  paragraph  presented  the  analogues  of  finite  element  h-  and 
p-refinerrient:  knot  insertion  and  order  elevation.  To  be  identical  to  h- 
refinement,  knot  insertion  must  be  performed  such  that  each  new  knot  has  a 
multiplicity  equal  to  the  polynomial  order  of  the  basis,  ensuring  C°  continuity 
at  each  knot.  Similarly,  if  we  begin  with  a  mesh  where  all  functions  are  already 
C°,  order  elevation  coincides  exactly  with  the  usual  notion  of  p- refinement 
in  the  finite  element  literature.  However,  knot  insertion  and  order  elevation 
provide  us  with  additional  possibilities  (see  [IT]). 

A  higher-order  alternative  elevation  strategy  present  itself;  it  takes  advantage 
of  the  fact  that  knot  insertion  and  order  elevation  do  not  commute.  If  a  unique 
knot  value  £  is  inserted  between  two  distinct  knot  values  in  a  curve  of  order 
p,  the  number  of  continuous  derivatives  of  the  basis  functions  at  £  is  p  —  1. 
If  we  subsequently  elevate  to  order  q,  the  multiplicity  of  every  knot  value  is 
increased  to  maintain  Cp_1  continuity  of  the  basis  functions  at  £.  If  instead  we 
order  elevate  the  original  curve  to  q  and  then  insert  a  unique  knot  value  £,  the 
basis  would  have  q  —  1  continuous  derivatives  at  £.  This  strategy  is  referred  to 
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as  k-refinement ;  it  has  no  analogue  in  standard  finite  element  analysis. 


Seven  C11  quadratic  functions 


knot  insertion 

Two  C°  linear  functions  ?  =  1/3. 2/3 


P  =  1  — »  p  =  2 


Fig.  5.  Non  commutativity  of  h-  and  p-refinement.  The  starting  point  is  one  linear 
element.  Classical  p-refinement  strategy:  knot  insertion  followed  by  order  elevation 
results  in  seven  C°  quadratic  functions.  New  higher-order  fc-refinement  strategy: 
order  elevation  followed  by  knot  insertion  results  in  five  C1  quadratic  functions. 


The  concept  of  ^-refinement  is  important  because  isogeometric  analysis  is  fun¬ 
damentally  a  higher-order  approach.  In  traditional  p-refinement  there  is  a  very 
inhomogeneous  structure  to  arrays  due  to  the  different  basis  functions  associ¬ 
ated  with  surface,  edge,  vertex,  and  interior  nodes.  Furthermore,  maintaining 
C°  continuity  during  the  refinement  process  implies  a  proliferation  in  the  num¬ 
ber  of  nodes,  In  /c- refinement,  there  is  a  homogeneous  structure  within  patches 
and  growth  in  the  number  of  control  variables  is  limited  (see  [17]). 


An  example  of  the  ^-refinement  process  and  its  comparison  with  traditional  p- 
refinement  is  given  in  Figure  5.  This  example  shows  that  fc-refinement  results  in 
fewer  functions  with  higher  continuity,  thus  fewer  control  variables  or  degrees 
of  freedom.  Starting  with  p  +  1  functions,  inserting  n  —  (p  +  1)  knots  (he.,  to 
obtain  n  basis  functions),  followed  by  r  order  elevations,  results  in  (r  +  l)n  —  rp 
Cp_1  basis  functions.  Using  fc-refinement,  that  is,  starting  with  p  +  1  functions, 
performing  r  order  elevations,  followed  by  the  insertion  of  n  —  (p  +  1)  knots, 
results  in  n  +  r  Cr+P~1  basis  functions.  More  details  on  this  can  be  found  in 
[17,  18,  32],  Keeping  in  mind  that  the  preceding  numbers  are  raised  to  the  d 
power  in  d  dimensions,  it  is  clear  that  the  fc-refinement  strategy  produces  fewer 
unknowns  than  p-refinement,  for  the  same  mesh  and  order  of  approximation. 


In  the  ^-refinement  process  on  a  fixed  mesh,  knots  are  added  at  the  boundary 
points,  increasing  their  multiplicity,  but  no  interior  knots  are  added.  It  is  inter¬ 
esting  to  note  that  in  the  periodic  case,  in  which  there  are  no  boundaries  and 
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consequently  no  open  knot  vectors  (see,  e.g.  Bazilevs  et  al.  [5]),  no  additional 
equations  are  engendered  by  ^-refinement. 


2.1.6  B-Spline  surfaces  and  solids 

Given  a  control  net  i  =  1  =  1  m  and  knot  vectors  5  = 

{£i, £n+p+i},  H  =  {rji, ...,  r)m+q+i},  a  tensor  product  B-Spline  surface  is 
defined  by: 


S(Z,  V)  =  E  E  (8) 

i=  1  j= 1 

where  NitP  and  Mhq  are  B-Spline  basis  functions  of  order  p  and  q  respectively. 
For  purposes  of  numerical  integration  of  arrays  constructed  from  B-Splines, 
“elements”  are  taken  to  be  knot  spans,  namely  [£* , C,;+  i ]  x  [VjiVj+ih  as  seen 
previously  in  the  one-dimensionnal  case  in  Figure  3.  As  presented  in  [32], 
integrals  are  pulled  back  to  the  parent  element  by  the  classical  change  of 
variables  formula  and  standard  Gaussian  quadrature  rules  are  employed  as 
in  finite  elements  (see  for  example  [30],  Chapter  3).  However,  we  note  that 
optimally  efficient  quadrature  rules  for  B-Spline  patches  do  not  seem  to  be 
known  yet. 

Analogously,  tensor  product  B-Spline  solids  can  be  defined;  given  a  con¬ 
trol  net  {Btjk},  i  =  1,.,,,  n,j  =  1, m,  k  =  1 ,...,/  and  knot  vectors 
‘  {si?  •••?  'Cn+p+l} ?  T~(-  {hi?  •••?  Vm+q+ 1}  Slid  Z  {Cl?  •••?  Ci+r+l}- 

n  m  l 

s&v,  0  =  EEE  NiiP(£)Mjtq(r})Lktr(()Bitjtk,  (9) 

i=l  j= 1  k= 1 

where  Ni>p,  Mj  q  and  Lk  r  are  B-Spline  basis  functions  of  order  p,q  and  r  re¬ 
spectively. 


2.2  Rational  B-Splines 


As  described  previously,  NURBS  are  built  from  B-Splines.  Specifically, 
NURBS  entities  in  can  be  built  from  a  projective  transformation  of 
B-Spline  entities  in  Md+1.  In  particular,  conic  sections,  such  as  circles  and 
ellipses,  can  be  exactly  constructed  by  projective  transformations  of  piecewise 
rational  quadratic  curves. 

To  obtain  a  NURBS  curve  in  lU,  we  start  from  Bf ,  i  =  1, ...,  n,  a  set  of  control 


12 


points  (“projective  points”)  for  a  B-Spline  curve  in  Md+1  with  knot  vector  5. 
The  control  points  for  the  NURBS  curve  are  defined  by: 


(Bi)j 


(B 


w\ 
i  >3 


Wj 


j  1, d, 


(10) 


where  ( Bt)j  is  the  jth  component  of  the  vector  Bi  and  wy  =  (Bf)d+ 1  is  referred 
to  as  the  ith  weight.  The  rational  basis  functions  and  NURBS  curve  are  given 
by: 


m) 


£?=1  NiA£)Wi' 


C(0  =  Y.Ri(0Bi. 

i=  1 


(11) 

(12) 


Rational  surfaces  and  solids  are  dehned  analogously  in  terms  of  the  rational 
basis  functions;  see  [32,  48,  52]  for  further  details: 


RTJ&V)  = 


Ni,P(0MjMwi,j 


£r=i  TZi  NiA0Mj,q(v)wij  ’ 


0  = 


■^i,p  )  Mj,q  ( Jl  )Rk,r  (C)  ^i,j,k 


E?=i  £7=i  £U  NiAtWjMLkAOwij* 


(13) 

(14) 


In  the  following,  we  summarize  noteworthy  properties  of  NURBS: 

•  NURBS  basis  functions  formed  from  an  open  knot  vector  constitute  a  par¬ 
tition  of  unity. 

•  The  continuity  and  support  of  NURBS  basis  functions  are  the  same  as  for 
B-Splines. 

•  NURBS  possess  the  property  of  affine  covariance. 

•  If  all  weights  are  equal,  NURBS  become  B-Splines. 

•  NURBS  surfaces  and  solids  are  the  projective  transformations  of  tensor 
product  piecewise  polynomial  entities. 


2.3  NURBS  as  a  basis  for  Isogeometric  Analysis 


The  isogeometric  analysis  framework,  based  on  NURBS  consists  of  the  follow¬ 
ing  features: 
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Finite  Element  Analysis 

Isogeometric  Analysis 

Nodal  points 

Control  points 

Nodal  variables 

Control  variables 

Mesh 

Knots 

Basis  interpolates  nodal 
points  and  variables 

Basis  does  not  interpolate 
control  points  and  variables 

Approximate  geometry 

“Exact”  geometry 

Polynomial  basis 

NURBS  basis 

Gibbs  phenomena 

Variation  dimishing 

Subdomains 

Patches 

Compact  support 

Partition  of  Unity 

Isoparametric  concept 

Affine  covariance 

Satisfy  patch  tests 

Table  1 

Comparison  of  Finite  Element  and  Isogeometric  Analysis  properties. 


•  A  mesh  for  a  NURBS  patch  is  defined  by  the  product  of  knot  vectors.  For 
example  in  three  dimensions,  a  mesh  is  given  by  S  x  7 -Lx  Z. 

•  Knot  spans  subdivide  the  domain  into  “elements.” 

•  The  support  of  each  basis  function  consists  of  a  small  number  of  elements. 

•  The  control  points  associated  with  the  basis  functions  define  the  geometry. 

•  The  isoparametric  concept  is  invoked.  The  coefficients  of  the  basis  functions 
are  the  degrees  of  freedom  or  control  variables. 

•  Three  different  mesh  refinement  strategies  are  possible:  analogues  of  h-  and 
p-refinement  and  the  new  higher-order  /e-refinement  scheme. 

•  Element  arrays  constructed  from  isoparametric  NURBS  are  assembled  into 
global  arrays  in  the  same  way  as  finite  elements  (see  [30],  Chapter  2). 

•  Dirichlet  boundary  conditions  are  applied  to  the  control  variables.  Homoge¬ 
neous  conditions  are  satisfied  pointwise.  For  inhomogeneous  conditions,  the 
boundary  values  must  be  approximated  by  functions  lying  in  the  NURBS 
space;  this  results  in  “strong”  but  approximate  satisfaction.  Another  option 
is  to  impose  Dirichlet  conditions  “weakly”  (see,  for  example,  [8]).  Neumann 
boundary  conditions  are  satisfied  naturally.  The  situation  is  very  similar  to 
standard  finite  elements  (see  [30],  Chapters  1  and  2). 
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Table  1  summarizes  the  similar  and  different  properties  of  isogeometric  anal¬ 
ysis  and  finite  element  analysis. 


3  B  formulation  for  linear  elasticity  using  NURBS 


3.1  The  boundary  value  problem  of  compressible  and  incompressible  linear 
elasticity 


The  boundary  value  problem  of  compressible  elasticity  for  a  body  fl  is  given 
by  the  following: 


Given  f  :  Q  — >  M3,  g  :  T9  — >•  M3,  and  h 

diver  +  f 

u 

ct  •  n 


:  r/(,  - 

->  M3,  find  u  :  Q  - 

-  R3  such  that: 

=  0 

in  H, 

(15) 

=  g 

on  Tg, 

(16) 

=  h 

on  Th, 

(17) 

n  is  the  exterior  unit  normal  on  T,  the  boundary  of  0,  g  is  the  prescribed 
displacement  on  T9  and  h  is  the  prescribed  traction  onT/,,  which  form  together 
the  boundary  T  =  T^  U  T9  of  ff,  and  f  is  the  body  force.  The  stress  tensor  cr 
is  defined  in  terms  of  the  strain  tensor  e  by  the  generalized  Hooke’s  law: 

„  1  „  T.  If  duj  du,  A 

£  —  V  u  —  ^(Vu  +  Vu  )o^«=2(^T  +  ^:J,  (18) 

(T  C  .  S  Or  (Jjj  Cijkl^kl  ( 1  ) 

The  Einstein  summation  convention  is  employed  for  spatial  indices  through¬ 
out. 

In  the  compressible  isotropic  linear  elastic  case,  Hooke’s  law  can  be  expressed 
in  terms  of  the  Lame  parameters  A  and  /j  by: 


Cijki  A<5 ijSki  T  p  (bf.kSji  T  biibji-'j , 

Up'  A  U  k ,  k  Oij  T  2  flc  j  j , 


where 


/z2/i  E 

(1  -  2u)  ’  ^  _  2(1  +  i/)’ 


and  v  is  Poisson’s  ratio  and  E  is  Young’s  modulus. 


(20) 

(21) 

(22) 


Clearly,  as  v  — »  |,  A  approaches  infinity.  The  value  v  —  \  thus  represents 
incompressibility.  The  constitutive  equation  needs  to  be  modified  in  this  case: 

Gij  pSij  Y  ( 2 3 j 
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where  p,  the  hydrostatic  pressure,  is  determined  as  part  of  the  solution  of 
the  boundary  value  problem.  As  p  represents  an  additional  unknown,  the 
kinematic  condition  of  incompressibility  must  be  introduced  as  an  additional 
equation: 

divu  =  Uk,k  =  0  in  H.  (24) 

However,  to  represent  the  nearly  incompressible  case,  that  is,  the  ratio  A//i 
is  large  but  not  infinite,  the  compressible  theory  can  be  applied  with  some 
modifications  in  the  discrete  case. 


3.2  Selective  and  reduced  integration 


Among  the  original  techniques  to  successfully  solve  the  nearly  incompressible 
case  are  selective  and  reduced  integration  techniques  as  presented  in  [28-30]. 
The  basic  idea  of  selective  integration  is  to  split  the  bilinear  form  of  the  varia¬ 
tional  formulation  into  its  A  and  p  contributions.  Using  vector  representation 
of  tensors,  and  introducing  the  strain- displacement  matrix  B  (see  [30] ,  chapter 
2  for  a  definition  of  B),  the  element  stiffness  matrix  is: 

ke  =  /  Br  DBdfl  (25) 

J  f2e 


The  material  properties  matrix  D  can  be  split  into  its  A  and  p  parts: 


D  =  //D  +  AD, 


(26) 


which  results  in  the  same  split  for  the  element  stiffness  matrix: 

ke  =  pke  +  Ake. 


(27) 


Since  X/p  3>  1,  the  numerical  values  of  the  entries  in  the  second  matrix  in 
Eq.  (27)  tend  to  be  very  large  compared  with  those  of  the  first.  The  idea  is  to 
use  a  standard  Gauss  quadrature  rule  to  compute  ke  and  a  reduced  (i.e.,  lower 
order  Gauss  quadrature)  rule  on  ke  to  reduce  the  number  of  incompressibility 
constraints.  The  simplest  example  is  the  piecewise  bilinear  quadrilateral  cle¬ 
ment  in  plane  strain  for  which  the  normal  rule  is  the  2x2  Gauss  points  rule 
and  the  reduced  rule  is  the  one-point  Gauss  rule.  Equivalence  theorems  be¬ 
tween  reduced/selective  integration  elements  and  mixed  formulation  elements 
have  been  obtained  by  Hughes  and  Malkus  [28,  35,  45].  These  theorems  show 
that  selective  and  reduced  integration  techniques  are  a  simple  way  of  attaining 
the  performance  of  mixed  methods  without  engendering  unwanted  degrees  of 
freedom.  This  seems  to  be  a  primary  reason  why  these  techniques  are  used  in 
large-scale  and  commercial  codes  [43,  44], 
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3.3  Strain  projection  and  the  B  method 


Some  drawbacks  of  the  selective  and  reduced  integration  rules  are  that  the 
equivalence  theorems  with  mixed  formulations  are  not  valid  in  the  axisym- 
metric  case  and  that  these  rules  are  difficult  to  generalize  in  anisotropic  cases. 
This  can  become  a  particular  problem  with  nonlinear  problems  since  the  tan¬ 
gent  moduli  always  exhibits  anisotropic  character.  These  difficulties  were  over¬ 
come  by  introducing  a  strain  projection  approach,  referred  to  as  the  B  method 
(Hughes  [29]). 

The  main  idea  in  the  strain  projection  approach  is  to  additively  split  the  strain 
tensor  into  its  deviatoric  and  dilatational  (he.,  volumetric)  parts 

e(u)  =  edev(u)  +  edil(u),  (28) 

where 

«*"(»)  =  l  (divu)  I  or  4“(u)  =  (29) 

and  I  is  the  identity  tensor. 

To  achieve  an  effective  formulation  in  the  nearly  incompressible  case,  the  di¬ 
latational  part  is  replaced  by  an  “improved”  dilatational  contribution  (e.g.,  a 
“projected”  one),  using  a  linear  projection  operator  7 r 

edll(u)  =  *  (edll(u))  .  (30) 


In  terms  of  B,  the  strain-displacement  matrix,  we  have 


with 


B  =  Bdev  +  Bdil, 

gdev  _  g  _  Bdil 


(31) 

(32) 


3.3.1  Weak  form  of  the  B  method 

Let  us  recall  the  usual  principle  of  minimum  potential  energy  to  obtain  the 
weak  form  in  linear  elasticity.  We  start  by  defining  the  trial  and  weighting 
spaces  S  =  {u  |  u  e  //1(H),u|rg  =  g}  and  V  =  {w  |  w  €  Lf1(h2),  w|r9  =  0}. 
Given  u  e  S  the  potential  energy  is  defined  by: 

n(u)  =  /  T(e(u  ))dkl  —  f  uf  dQ  —  f  u  •  hdT,  (33) 

Jo.  Jn  Jrh 
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where  T  is  the  strain  energy  density: 


^t(e(u))  =  -CijfliSij(u)£ki(u). 


(34) 


Minimizing  the  potential  energy  is  equivalent  to  satisfying  the  variational 
equation,  or  weak  form;  that  is  given  u  G  S,  w  €  V  and  eel, 


where 


ddhu  +  ew) 
de 


„o=  o  «  o(w,  u)  =  (w,  f )  +  (w,  h)r,, 


(35) 


a(w,  u)  = 

/  £ij(w )cijki£ki(u)dQ, 

(36) 

Jn 

/> 

(w,f)  = 

/  u  •  f dfl, 

(37) 

Jn 

/> 

w>h)rh  = 

/  u  •  hdT. 

Jrh 

(38) 

Considering  now  the  potential  energy  in  terms  of  the  modified  strain  tensor 
e: 

fl(u)  =  /  ^(e(u ))d£l  —  I  u  •  fe/fl  —  f  u  •  hdT, 

Jn  Jn  Jrh 

and  recalling  that  e  is  a  linear  operator,  one  can  write: 

e(u  +  ew))  =  e(u)  +  ee(w). 


(39) 


(40) 


Therefore,  the  directional  derivative  of  the  modified  strain  energy  density  is 
given  by: 

<9T(£(u  +  ew)) 


de 


E=0- 


de. 


-(e(u))  £ij{ w), 


which  allows  us  to  define  the  modified  stress  tensor  as: 

&ij  CijkiEki. 

The  minimization  of  the  modified  potential  energy  finally  becomes: 


(41) 


(42) 


0  =  /  £ij{w)cijki£ki(u)dn  -  /  W  •  fciO  -  /  w  ■  hdT,  (43) 

Jn  Jn  Jrh 

which  represents  the  variational  equation  for  the  B  method.  The  first  integral 
in  Eq.  (43)  is  a  bilinear  form,  and  expressing  it  the  following  way, 


a(w,u)=  /  £ij(w)cijki£ki(u)dn, 
Jn 


(44) 


18 


the  B  variational  formulation  can  be  stated  as: 

Find  uGiS,  such  that  Vw  €  V 
8(w,u)  =  (w,f)  +  (w,h)r„. 


(45) 


3-4  B  and  higher-order  NURBS 

3. 4-1  Projection  operator  and  space  with  NURBS 

The  use  of  the  B  method  within  isogeometric  analysis  requires  further  inves¬ 
tigation  into  the  choices  of  the  projection  operator  and  the  associated  space 
onto  which  the  projection  will  be  performed.  Since  the  technique  has  been 
applied  mostly  to  piecewise  bilinear  and  trilinear  finite  elements,  and  we  want 
to  make  intensive  use  of  the  properties  of  high-order  fc-refined  NURBS,  these 
topics  need  to  be  studied  without  any  assumption  on  the  order  of  approxima¬ 
tion. 

In  the  discrete  case,  we  have 

n 

Uh(x)  =  ^uAlVA(x),  (46) 

A= 1 

likewise 

n 

Wft(x)  =  ^wAlVA(x),  (47) 

A= 1 

where  NA  are  the  NURBS  basis  functions  and  u"4  and  wA  are  the  associated 
control  variables.  This  reduce  to  the  usual  shape  functions  and  degrees  of 
freedom  when  classical  finite  elements  are  utilized. 

In  developing  the  B  method  for  higher-order  finite  elements  and  NURBS, 
we  need  to  define  the  linear  projection  operator  and  the  spaces  upon  which 
to  project  the  dilatational  strain.  Throughout,  we  use  the  L2  projection  of 
the  strains.  For  the  spaces,  we  define  the  following  procedure:  assume  the 
displacement  space  is  given.  We  shall  refer  to  it  as  Qp,  that  is  quadrilateral,  or 
hexahedral,  elements  of  order  p.  The  continuity  of  Qp  elements  within  a  patch 
can  be  any  order  k  from  0  to  p—  1.  In  this  paper,  we  are  particularly  interested 
in  elements  of  maximal  continuity,  namely,  Ck  =  Cp_1.  We  always  assume  an 
open  knot  vector  construction  so  that  only  C°  continuity  is  attained  across 
patch  interface.  The  basis  functions  for  the  projected  dilatational  strain  are 
taken  to  be  one  order  lower,  and  usually  one  order  of  continuity  lower,  namely 
the  space  Qp-i,  of  continuity  Ck  =  Cp~2.  The  only  exception  occurs  when 
p  >  2,  but  there  are  lines  or  surfaces  of  C°  continuity  within  a  patch.  In  this 
situation,  the  projected  space  is  also  taken  to  be  C°  continuous  across  those 
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lines  or  surfaces.  There  is  nothing  fundamental  about  this  choice.  It  is  simply 
a  convenience  due  to  the  data  structure  we  employ  in  our  code.  In  practice, 
C°  lines  or  surfaces  within  a  patch  exist  when  conic  sections  such  as  circles, 
are  constructed  by  standard  NURBS  algorithms  (see  e.g.,  Piegl  and  Tiller  [48] 
and  Hughes  et  al.  [32]). 


As  an  example  of  what  our  construction  produces,  consider  the  displacement 
space  Q\ .  This  is  the  space  of  bilinear  quadrilaterals,  or  trilinear  hexahedra, 
and  is  C°  continuous  across  element  boundaries  (which  correspond  to  knots  in 
this  case).  The  space  for  projected  dilatational  strain  is  then  Q0,  of  continu¬ 
ity  class  C~l,  that  is,  piecewise  constants.  This  element  becomes  the  classical 
mean  dilatational  element  (see  Hughes  and  Allik  [31],  Nagtegaal  et  al.  [47] 
and  Hughes  [30])  referred  to,  herein,  as  Qi/Qo-  We  wish  to  emphasize  that  we 
recognize  that  there  is  no  guarantee  that  our  construction  produces  discrete 
approximations  satisfying  the  LBB  condition  (see  e.g.,  [10,  11,  62]).  In  fact, 
it  is  entirely  obvious  that  this  is  the  case,  the  Qi/Qo  element  being  a  prime 
example  of  an  element  that  fails.  However,  despite  the  theoretical  deficiency, 
it  must  be  noted  that  Qi/Qo  is  without  doubt  the  most  utilized  element  in  the 
history  of  nonlinear  structural  finite  element  analysis,  and  it  is  widely  used 
in  fluids  as  well.  There  are  other  practical  issues  at  play  here  (see  Hughes 
[30]).  We  will  say,  however,  that  our  construction  produces  a  good  balance 
between  the  number  of  displacement  degrees  of  freedom  and  number  of  di¬ 
latational  constraints  on  each  patch.  Asymptotically,  “constraint  ratios”  are 
2  in  two  dimensions  and  3  in  three  dimensions,  in  all  cases  (see  Hughes  [30] 
for  a  discussion  of  constraint  ratios).  But  from  the  LBB  perspective,  we  may 
have  slightly  too  many  constraints,  at  least  for  some  cases.  We  note  that  for 
the  discretizations  of  main  interest  herein,  namely  p  >  2,  Qp/Qp_i,  of  inter¬ 
nal  patch  continuity  Cp_1/Cp_2,  there  are  no  known  theoretical  results,  one 
way  or  the  other,  at  this  point  in  time.  However,  results  from  the  p-method 
community  seem  to  indicate  that  incompressible  locking  is  not  an  issue  for  suf¬ 
ficiently  great  p  (see  Suri  [61],  Szabo  and  Babuska  [62]  and  Szabo  et  al.  [63]). 
The  construction  we  advocate  takes  a  middle  road.  It  does  not  reduce  the 
number  of  constraints  sufficiently  to  satisfy  the  LBB  condition  for  all  cases, 
but  it  does  reduce  it  to  a  value  that  seems  to  attain  a  good  balance  in  all 
cases.  Our  experience  with  it  is  good,  as  the  numerical  results  will  indicate 
in  the  sequel.  On  the  other  hand,  our  experience  with  higher-order  elements 
without  projection  is  that  they  are  overconstrained  and  not  competitive  with 
the  reduced  constraint,  projection  elements  advocated.  Likewise,  spectral  ele¬ 
ments  of  order  Qp/Qp_2,  and  continuity  class  C° fC~x  are  known  to  satisfy  the 
LBB  condition  [39],  but  convergence  is  suboptimal  and  the  minimum-order 
element  is  In  fairness,  the  spectral  element  approach  is  focused  on 

very  high-order  elements  and  high  precision  applications  exclusively,  not  the 
low-order  end  of  the  spectrum.  Nevertheless,  in  nonlinear  structural  analysis, 
low-order  elements  are  extremely  important. 
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Qp 


Qp- 1 


Fig.  6.  Basis  functions  Qp/Qp- 1,  p  =  1,2,3,  for  a  one-dimensional  patch  of  four 
elements.  All  cases  attain  continuity  Cp_1 /Cp~2 . 

We  will  now  give  the  technical  details  of  the  construction  of  the  spaces.  On 
each  patch,  open  knot  vectors  are  employed.  Tensor  product  constructs  are 
utilized  so  we  focus  on  the  situation  in  each  direction  separately.  The  B-Spline 
basis  is  completely  determined  from  Eqs.  (2)  and  (3).  We  need  to  specify  the 
order  of  the  space  and  the  knot  vector.  We  begin  with  the  displacement  space, 
assumed  to  be  of  order  p  >  1.  The  knot  vector,  denoted  Sp,  is  assumed  to 
have  the  following  form, 


=  {(MW),  Eint,  MW},  (48) 

p-\- 1  copies  p+1  copies 

where,  for  simplicity,  we  have  assumed  the  initial  and  final  knots  are  located 
at  0  and  1,  respectively.  The  multiplicities  of  the  initial  and  final  knots  are 
p  +  1.  Eint  denotes  the  vector  of  internal  knots.  Each  internal  knot  may  have 
a  different  multiplicity,  allowable  values  being  1,2,3,  ...,p.  The  case  we  are 
primarily  concerned  with  in  this  paper  is  each  internal  knot  having  multiplicity 
1  which  results  in  maximal  smoothness  of  continuity  class  Cp~l  on  each  patch. 
The  corresponding  knot  vector  for  the  projected  space,  denoted  Ep_i,  is  given 
by 

W  =  {MW,  Eint,  MW}-  (49) 

p  copies  p  copies 

The  order  of  the  projected  space  is  taken  to  be  p  —  1  >  0.  Note  that  the 
internal  knot  vector  is  exactly  the  same  as  for  the  displacement  space,  whereas 
the  initial  and  final  knots  have  multiplicity  p.  The  span  of  the  projected  space 
is  precisely  the  span  of  the  derivatives  of  all  functions  in  the  displacement 
space.  An  example  of  the  spaces  Qp/Qp- 1,  p  —  1,  2,  3,  in  the  general  case,  for 
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a  one-dimensional  patch  of  four  elements,  is  given  in  Figure  6.  We  see  that 
Cp~l /Cp~2  continuity  is  achieved  in  all  cases. 


Qp  Qp—i 


Fig.  7.  Basis  functions  Qp/Qp~  1,  p  =  2,3,  with  a  point  of  C°  continuity  within 
the  patch,  for  a  one- dimensional  patch  of  four  elements.  All  cases  attain  continuity 
Cp_1/Cp_2,  except  at  the  repeated  knot  £  =  0.75  where  the  continuity  is  only  C°. 

The  exception  to  the  general  case  occurs  when  there  are  lines  or  surfaces  of  only 
C°  continuity  within  a  patch,  as  mentioned  previously.  Let  us  assume  that 
has  one  or  more  knots  having  multiplicity  p,  signifying  C°  continuity.  Then,  in 
the  space  Sp_!,  Eint  needs  to  be  replaced  with  Emt ,  which  is  identical  to  Eint 
except  for  the  knots  having  multiplicity  p;  in  Eint  these  knots  have  multiplicity 
p  —  1,  preserving  C°  continuity  of  the  projected  space  within  each  patch.  An 
example  of  the  spaces  Qp/Qp- 1,  p  =  2,  3,  with  a  point  of  C°  continuity  inside 
the  patch,  is  given  in  Figure  7.  Only  C°  continuity  is  achieved  across  the 
repeated  knot  £  =  0.75,  and  Cp_1/Cp_2  continuity  is  achieved  elsewhere  on  the 
patch  interior. 


The  construction  of  higher-dimensional  B-Splines  proceeds  by  assigning  the 
polynomial  order  and  knot  vector  for  each  dimension.  For  example,  in  three 
dimension,  we  have,  x  Hq  x  Zr  and  Sp_i  x  Hq-\  x  Zr_ i,  where 


nq  =  {0, 0, ...,  O.Hint,  1, 1, ...,  1}, 

J  V  V  V  V 

^  ^  N - V - ' 

q+ 1  copies  q+ 1  copies 

nq.x  =  { 010^20 ,  Hint ,  , 

q  copies  q  copies 

Zr  =  {0, 0, ...,  0,  ZinU  1, 1, ...,  1}, 


’V' 

r+1  copies 


■  v*" 

r+1  copies 


Zr~ i  —  (0,  o -j •  • ,  0,  Zjnt,l,  1, ...,  1}, 

r  copies 


v 

r  copies 


(50) 

(51) 

(52) 

(53) 


and  the  internal  knots  in  Tiint  and  Zint  may  take  multiplicities  in  the  range 
1,  2,  3, ...,  q ,  and  1,  2,  3, ...,  r,  respectively.  As  before,  if  there  are  knots  in  7 iint 
and  Zint  having  multiplicity  q  and  r,  then  'Hmt  and  Zint  need  to  be  replaced 
by  'Hint  and  Zmt ,  constructed  in  a  similar  fashion  as  to  Sint. 
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In  practice,  the  NA,s  are  constructed  as 

NA(x.)  =  Na  o  0-1(x),  (54) 

where  {N}  is  a  NURBS  basis  in  parametric  coordinates  and  (f)  is  the  exact 
mapping  between  the  parametric  and  geometric  domains.  We  construct  the 
“tilde  basis”  {N},  which  corresponds  to  the  projection  space,  by 

Na  (x)  =  Na°  0_1(x),  (55) 

where  {N}  is  the  lower  order  NURBS  basis  built  on  the  same  parametric 

domain.  As  described  previously,  we  take  the  NA’s  to  be  one  order  lower  than 
the  Na,s  to  reduce  the  number  of  incompressibility  constraints. 

Note  that,  even  in  the  case  of  lowest-order  elements  (he.,  bilinear  and  trilin- 
ear),  we  still  use  the  exact  geometrical  mapping.  This  means  our  lowest-order 
elements  are  isogeometric  and  precisely  fit  curved  boundaries.  We  believe  that 
Barth  [3]  was  the  first  to  use  this  approach,  and  to  demonstrate  its  effective¬ 
ness  in  compressible  fluid  calculations. 

In  the  discrete  case,  Eq.  (30)  becomes: 

n 

4‘V)  =  £v<4  (se) 

A=1 


where 


eA  = 
ij 


(nb,  egV))  =  /  NBef(uh)<m,  (57) 


B= 1 


B= 1 


that  is 


n  n  p 

“'*)  -  I  V  Vt M-/b  /  NB 

A,B=1C=1 


~dil  j 


dNc 

dxu 


dQ  u^dij, 


(58) 


and  M  is  the  “mass”  matrix  of  the  tilde  basis,  namely 


Mab=  Na,Nb  =  NANBdn. 


n 


(59) 


In  summary,  the  procedure  corresponds  to  L 2  projection  of  eff  onto  the  { A^} 
basis. 


3-4-2  Implementational  aspects 

We  consider  aspects  of  solving  the  global  matrix  system  in  this  section. 
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For  an  isotropic  homogeneous  linear  elastic  material,  with  Hooke’s  law  given 
by  Eq.  (21),  the  discrete  version  of  the  bilinear  form  of  Eq.  (44)  is  given  by: 


a(wh,uh)  = 


i>2  (wt  [  NjcijkiNfdQ  u% 

tsL i  V  Jn 


A,B= 

1 

+  3 


(3A  +  2k)  V  wf  (JVf,  Nc)MpD(ND,  JVf )  u 

C,D=  1 


B 
k  I  ) 


where 


&ijkl  (  &ik&jl  “I-  SU5jk  ^ijSkl 


(60) 


(61) 


Due  to  the  inverse  of  M,  the  second  term  in  Eq.  (60)  increases  the  popula¬ 
tion  of  the  stiffness  matrix  on  each  patch  for  p  >  2.  Note,  we  always  assume 
use  of  patches  constructed  with  open  knot  vectors.  This  means  that  the  dis¬ 
placement  field  is  continuous  across  patch  interfaces,  but  no  smoother.  Conse¬ 
quently,  the  tilde  basis  will  be  discontinuous  across  patches  and  M-1  will  be 
uncoupled  from  patch  to  patch.  Nevertheless,  if  we  use  a  direct  solver  to  solve 
the  global  equation  system,  we  need  to  account  for  increased  coupling  of  the 
equations  due  to  M-1  on  each  patch.  There  are  at  least  two  ways  to  circum¬ 
vent  the  effect  of  the  increased  coupling.  One  is  to  use  an  iterative  strategy 
that  does  not  require  the  assembly  of  the  stiffness  matrix,  such  as  Conju¬ 
gate  Gradients,  to  solve  the  global  problem.  Within  each  conjugate  gradient 
iteration,  a  direct  solver  can  be  used  to  evaluate  M_1  patch- wise,  retaining 
its  sparse  band-profile  structure.  This  procedure  can  be  used  to  solve  very 
large  problems.  We  have  used  it  extensively  in  our  calculations  and  found  it 
to  be  very  efficient.  A  second  possibility  is  to  replace  M  with  a  diagonal,  or 
“lumped”  approximation.  This  would  only  need  to  be  done  in  the  left-hand- 
side  matrix,  and  so  would  be  interpreted  as  a  preconditioner.  In  this  case,  the 
band-profile  structure  of  the  preconditioner  would  be  only  slightly  larger  than 
for  the  system  constructed  without  projection.  Using  the  consistent  M  on  the 
right-hand-side  would  ensure  the  full  accuracy  of  the  projection  procedure. 
Convergence  would  require  one  or  more  iterations,  but  this  involves  only  a 
forward  reduction  and  back  substitution  for  each  additional  iteration  with  an 
existing  factorized  array  when  employing  a  direct  solver. 


3.5  Numerical  Examples 


We  consider  NURBS  displacement /projection  spaces  Qp/Qp-i  with  p  = 
1,2, 3, 4.  All  internal  knots  have  multiplicity  1,  resulting  in  maximal  smooth¬ 
ness,  that  is  Cp~x/Cp~2.  For  fixed  p,  we  examine  “h- refinement,”  that  is,  we 
subdivide  the  mesh  by  knot  insertion.  For  a  fixed  mesh,  we  also  consider 
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^-refinement,  that  is,  we  increase  p  before  inserting  knots  to  obtain  the  desired 
mesh,  ^-refinement  produces  very  few  additional  degrees  of  freedom  (he., 
control  variables)  as  no  internal  knots  are  added  and  only  the  multiplicities 
of  the  initial  and  final  knots  are  increased  by  1.  This  results  in  one  layer  of 
additional  degrees  of  freedom  per  dimension.  For  example,  in  a  one-patch 
three-dimensional  mesh  of  3n3  degrees  of  freedom,  one  level  of  ^-refinement 
results  in  (3 n  +  3)3  degrees  of  freedom  (see  Cottrell  et  al.  [17]  for  further 
details). 

In  the  following  we  present  two  numerical  examples  for  nearly  incompress¬ 
ible  linear  elasticity  to  demonstrate  the  proposed  approach.  Both  are  plane 
strain  two-dimensional  problems  but  are  solved  using  a  three-dimensional  for¬ 
mulation  with  two  linear  functions  in  the  out-of-plane  direction  and  boundary 
conditions  enforcing  the  plane  strain  constraint. 


3.5.1  Cook’s  membrane 


E  =  240.565  MPa 
v  =  0.4999 
F  =  100  N/mm 


Fig.  8.  Geometry,  loading,  boundary  conditions,  material  parameters  and  quantity 
of  interest  for  the  plane  strain  Cook’s  membrane. 

This  problem  has  been  solved  by  many  authors  to  test  nearly  incompressible 
formulations  under  combined  bending  and  shear  (see,  e.g.,  [14,  16,  19,  40]).  A 
tapered  panel  is  clamped  on  one  side  and  subjected  to  a  uniform  shear  load 
on  the  opposite  side.  The  other  sides  are  traction  free.  The  geometry,  loading 
and  boundary  conditions  are  given  in  Figure  8.  The  quantity  of  interest  is 
the  vertical  displacement  of  the  top  right  corner  of  the  plate  as  presented  in 
Figure  8.  Its  convergence  as  a  function  of  the  number  of  elements  per  edge  is 
usually  used  as  a  criterion  in  investigating  performance.  The  results  are  plotted 
in  Figure  9  for  the  four  orders  of  approximation  with  the  standard  Galerkin 
approximation,  without  any  treatment  of  the  incompressibility  condition,  and 
with  the  proposed  B  method.  We  observe,  as  expected,  that  the  piecewise 
bilinear  quadrilateral  element  suffers  from  severe  locking.  Increasing  the  poly- 
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Fig.  9.  Cook’s  membrane.  Vertical  displacement  of  top  right  comer  versus  number 
of  elements  per  edge  with  and  without  B  for  various  NURBS  orders  obtained  from 
^-refinement. 


Number  of  Degrees  Of  Freedom 


Fig.  10.  Cook’s  membrane.  Vertical  displacement  of  top  right  corner  versus  number 
of  degrees  of  freedom  with  and  without  B  for  various  NURBS  orders  obtained  from 
^-refinement.  The  case  “No  B  Ql”  is  offscale. 
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nomial  order  of  the  NURBS  functions  improves  the  results,  however,  we  can 
see  that  C 1  quadratic  NURBS  still  suffer  from  locking  for  coarse  meshes.  The 
use  of  B  significantly  improves  the  results,  and  good  convergence  is  attained 
even  with  Qi/Qq.  Note  that  the  results  obtained  with  Qi/Qo  are  better  than 
Q2  and  that  Q2IQ1  gives  better  results  than  Q4.  Finally,  we  observe  that  with 
B  Q4/Q3  we  obtain  converged  results  with  only  a  mesh  of  four  elements. 

It  is  also  interesting  to  study  convergence  as  a  function  of  the  number  of  un¬ 
knowns.  The  results  are  shown  in  Figure  10,  where  a  logarithmic  scale  has 
been  used  for  the  abscissa.  Note  here,  that  in  terms  of  unknowns,  Q1/Q0  per¬ 
forms  almost  the  same  as  Qs  and  Q4,.  We  also  see  that  the  results  with  Q2/Q1 
and  Q3/Q2  are  virtually  equivalent  in  terms  of  unknowns  and  the  best  per¬ 
formance  is  attained  with  Q4/Q3.  These  residts  demonstrate  the  superiority 
of  high-order  /j-rehned  NURBS:  increasing  the  polynomial  order  and  the  con¬ 
tinuity  only  slightly  increases  the  number  of  unknowns,  which  enables  us  to 
obtain  very  good  results  with  only  a  few  C3  quartic  functions:  for  Q4/Q3,  4 
elements  and  36  degrees  of  freedom  is  enough  to  obtain  the  reference  value, 
whereas  we  need  at  least  1024  elements  and  1089  degrees  of  freedom  to  obtain 
about  the  same  results  with  Qi/Qq. 


3.5.2  Infinite  plate  with  circular  hole  under  in-plane  tension 

The  next  problem  is  a  plane-strain  infinite  plate  with  a  hole  under  tension. 
This  problem  has  been  studied  previously  with  isogeometric  analysis  in  [32] 
assuming  an  isotropic  compressible  linear  elastic  medium.  It  is  interesting 
from  the  geometrical  point  of  view  because  quadratic  NURBS  can  exactly 
represent  the  circular  hole  and  the  existence  of  an  analytical  solution  allows 
us  to  focus  on  the  convergence  rates  that  the  proposed  method  can  attain 
without  geometrical  approximation.  In  the  nearly  incompressible  regime,  this 
problem  has  been  of  less  interest  than  the  previous  one.  Some  results  in  the 
incompressible  limit  using  meshless  methods  can  be  found  in  [23,  27].  The 
infinite  plate  is  modeled  by  a  quarter  plate.  The  exact  solution  (Timoshenko 
and  Goodier  [64],  pp.  90-93)  is  evaluated  at  the  boundary  of  the  quarter  plate 
and  applied  as  a  Neumann  boundary  condition. 

The  exact  solution  is  given  by  the  following: 

MM)  =  y(!  -  ~j)  +  y(x  ~~  +  3~r)  cos2^>  (62) 

<roo(r,0)  =  y(l  +  y)  -  y(l  +  3y)cos2d,  (63) 

<7ro{r,0)  =  — y  (1  +  2y  -  3y)sin2 9.  (64) 

The  geometry,  loading,  boundary  conditions  and  parameters  are  shown  in 
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Fig.  11.  Geometry,  loading,  boundary  conditions  and  material  parameters  for  the 
infinite  plate  with  a  hole  under  tension. 

Figure  11.  R  is  the  radius  of  the  hole,  L  the  length  of  the  finite  quarter 
plate,  and  Tx  is  the  magnitude  of  the  applied  stress  for  the  infinite  case. 
The  value  for  v  was  chosen  to  be  very  close  to  the  incompressible  case  to 
study  the  convergence  of  the  B  isogeometric  analysis  in  that  case.  A  rational 
quadratic  basis  is  the  minimum  order  capable  of  representing  the  geometry. 
The  sequence  of  meshes  obtained  from  /i-refinement  used  for  the  convergence 
study  are  shown  in  Figure  12. 


Fig.  12.  Infinite  plate  with  a  hole:  sequence  of  meshes  produced  by  /i-refinenrent. 

The  geometry  of  the  coarsest  mesh  can  be  exactly  represented  by  quadratic 
NURBS.  It  consists  of  two  elements,  shown  on  the  left  of  Figure  12,  and  is 
defined  by  the  knot  vectors  S2  and  7i2,  from  which  the  knot  vectors  Si  and 
Hi  are  constructed  to  build  the  corresponding  linear  basis: 

Z2xH2  =  {0,  0, 0,  0.5, 1, 1, 1}  x  {0,  0,  0, 1, 1, 1},  (65) 

Sx  x  Hi  =  {0, 0, 0.5, 1, 1}  x  {0, 0, 1, 1}.  (66) 


We  define  the  relative  error  as  the  error  normalized  by  the  corresponding  value 
of  the  exact  solution.  Convergence  results  for  the  relative  error  in  the  L2  norm 
of  displacement,  relative  error  in  the  L2  norm  of  stress,  and  relative  error  in 


mesh  parameter  h 


Fig.  13.  Plate  with  a  hole.  Convergence  curves  in  relative  L2  norm  of  displacement, 
relative  L2  norm  of  stress  and  relative  energy  norm  with  and  without  B  for  various 
NURBS  order  obtained  from  /c-refinement. 
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the  energy  norm  are  shown  in  Figure  13.  The  cubic  and  quartic  NURBS  are 
obtained  from  ^-refinement  of  the  coarsest  quadratic  mesh.  The  mesh  param¬ 
eter  h  is  defined  as  the  maximum  distance,  in  the  physical  space,  between 
diagonally  opposite  knot  locations.  As  can  be  seen,  the  B  method  obtains 
good  convergence  rates  in  the  three  norms  with  relatively  coarse  meshes.  Note 
that  the  standard  displacement  based  formulation  performs  relatively  poorly, 
especially  in  the  L2  norm  of  stress,  and  needs  relatively  fine  meshes  to  attain 
convergence  rates  equivalent  to  what  is  obtained  with  the  B  formulation.  Even 
when  seemingly  optimal  asymptotic  rates  of  convergence  are  attained  in  the 
case  of  the  standard  Qp  elements,  the  error  is  orders  of  magnitude  greater 
than  for  the  corresponding  projected  Qp/Qp~i  B  elements.  For  example,  in 
stress  the  difference  is  four  orders  of  magnitude!  This  result  clearly  shows  that 
optimal  rate  of  convergence  is  not  the  only  issue  to  be  considered. 


4  F  formulation  for  nonlinear  elasticity 


Jhl  Constitutive  equations 


The  central  idea  in  the  development  that  we  propose  is,  as  done  in  the  geo¬ 
metrically  linear  case,  to  split  the  tensor  that  measures  the  deformation  into 
its  deviatoric  (volume  preserving)  and  volumetric-dilatational  parts.  In  the  fi¬ 
nite  deformation  case,  the  deformation  gradient  F  is  the  relevant  tensor,  and, 
contrary  to  the  small  deformation  case,  the  split  is  multiplicative  rather  than 
additive.  This  multiplicative  decomposition  has  been  exploited  previously  by 
Flory  [26],  Hughes  et  al.  [38],  Sirno  et  al.  [54,  60]  (within  a  three  held  Hu- 
Washizu  principle),  and  more  recently  by  de  Souza  Neto  et  al.  [20,  22]  in  an 
alternative  F  approach.  The  work  we  present  here  shares  features  with  these 
techniques.  We  want  it  to  be  a  simple  pure  displacement  formulation,  but 
having  greater  generality  and  a  more  rigorous  theoretical  base  than  previous 
approaches. 

We  start  by  writing  the  multiplicative  split  of  the  deformation  gradient: 

F  =  FdilFdcv,  (67) 

where 

detF  =  J  =  detFdil  and  detFdev  =  1,  (68) 

which  leads  to: 

Fdev  =  J-V3F  and  Fdil  =  J1/3I.  (69) 

We  now  define  a  modified  deformation  gradient  F  in  terms  of  the  deviatoric 
part  of  the  deformation  gradient  Fdev  and  a  modified  dilatational  part  of  the 
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deformation  gradient  F' 1  : 


(70) 


jp  _  jpdiljpdev 

where 

FdiI  =  vr(Fdil)  =  tt(  J1/3)I  =  JV3 1,  (71) 

with  7T  a  linear  projection  operator  identical  to  the  one  proposed  previously 
for  the  linear  case.  We  can  finally  write,  keeping  the  “bar”  notation  for  the 
projection: 


ccF. 

(72) 

J1/3 

(73) 

J1/3' 

We  assume  hyperelastic  homogeneous  material  behavior  for  which  there  exists 
a  free-energy  function3  T  that  depends  on  the  Cauchy-Green  tensor  C  =  F1  F 
and  from  which  the  second  Piola-Kirchhoff  stress  tensor  is  derived  as: 


S 


(0*(C) 
'  dc 


(74) 


The  standard  additive  decomposition  of  T  (see  for  example  [58])  into  a  volu¬ 
metric  part  depending  only  on  J  and  an  isochoric  part  is  used: 


'Id./.  C  j  —  'ldm!'./  i  +  C  ) 


(75) 


4-2  Variational  equation  and  weak  form 


Let  B  be  the  reference  configuration  of  a  body  and  B'  be  its  deformed  config¬ 
uration.  We  introduce  a  mapping  <fi  :  B  — >  B'  that  takes  a  point  X  G  B  to  a 
point  x  =  </>(X)  G  B' .  The  displacement  of  a  particle  from  its  initial  position 
X  to  its  current  position  x  is  given  by 

u(X)  =  </>(X)  -  X  =  x  -  X  (76) 

The  deformation  gradient  is  the  second  order  tensor  defined  by: 

F(X)  =  VA>(X)  =  (77) 


The  boundary  value  problem  for  finite  deformation  elasticity  for  a  body  with 


3  also  called  stored  energy  or  strain  energy  function 
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reference  configuration  B  is  given  by  the  following: 


Given  T  :  B  — >  M3,  <p  '■  d<t>B  —> ►  M3,  and  H  :  d^B  — >  M3,  find  <p  G  S  such  that: 

Vx  •  P  +  T  =  0  in  B,  (78) 

<p  =  cj)  on  d<pB,  (79) 

P  •  N  =  H  on  dnB ,  (80) 

where  N  is  the  unit  exterior  normal  on  d B,  the  boundary  of  B,  (ft  is  the 
prescribed  deformation  on  d^B  and  H.  is  the  prescribed  Piola  traction  on  d^B, 
which  form  together  the  boundary  d B  =  d^B  U  d^B.  T  is  the  body  force  per 
unit  initial  volume,  S  is  the  space  of  admissible  deformations  defined  by 

5  =  {4>  :  -*•  M3  |  det(VA^)  >  0,  <p\a<j>1&  =  0},  (81) 

and  P  is  the  first  Piola-Kirchhoff  stress  tensor: 

P  =  FS.  (82) 


Corresponding  to  what  was  done  previously  in  the  geometrically  linear  case, 
given  the  mapping  (p  G  S,  the  potential  energy  is  given  by 

11(0)  =  [  T(E (0))  dV  -  [  <p  ■  T  dV  -  [  p-Hdr,  (83) 

Jb  Jb  J  dnB 

where  E  =  ^(F7  F  — I)  is  the  Green- Lagrange  strain  tensor.  The  stationarity  of 
the  potential  energy  leads  to  the  variational  equation  for  the  finite  deformation 
case  (see  for  example  [58]). 

We  now  introduce  the  modified  potential  energy 


fl(0)  =  /  (E(0))  dV  —  [  <p-T  dV  -  [  p-Hdr,  (84) 

Jb  Jb  JdnB 

where  the  modified  Green-Lagrange  strain  tensor  is  defined  in  terms  of  the 
modified  deformation  gradient: 

E  =  I(FtF-I).  (85) 

We  write  the  stationarity  of  the  modified  potential  energy: 


011(0 +  eW) 

de 


6  =  0  — 


0, 


where  <p  E  S,  W  G  V  =  {W  :  B  — ►  M3  |  W  =  0}  and  eel. 


(86) 
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We  need  to  compute  the  directional  derivative  of  the  strain  energy  density, 
keeping  in  mind  that  the  projection  operator  is  linear: 


<9T(E(<7>  +  eW)) 
We 


h=o= 


0tf(E) 


dE 


ij 


(E(0)) 


dEu((f)  +  e  W) 
de 


h=o 


(87) 


The  calculation  of  the  directional  derivative  of  the  modified  Green- Lagrange 
strain  tensor  is  presented  in  Appendix  A,  the  result  is: 


dE(<j b  +  eW) 
We 


E=0  =  DEu{<f>)  ■  W 

1  fn(J1/3Vx  •  W) 
JIT! 

(FtVaW)  , 


Vx  •  W  F  F 


“T  O' 

where  a  is  dehned  in  Eq.  (73). 

We  then  define  the  modified  second  Piola-Kirchhoff  stress  tensor 

s,m  =  s„(  e(0»  =  7^(E(<w). 


(88) 


d  E 


(89) 


ij 


and  Eq.  (87)  becomes 


SnDEutf)  ■W  =  SJ 


ij 


1  /^(J^V*  •  W) 

7177 


Vs  ■  W  )  FJ'F  +  a(F'J  VAW) 

(90) 


where  Va’  •  W  =  tr(VAWF  x)  (see  Eq.  (B.15)  in  Appendix  B). 
Defining  the  modified  first  Piola-Kirchhoff  stress  tensor  as 

Pu  =  FijSu , 

and  introducing  the  following  notation 


Wij  =  aWij  +  - 


1  ■  W) 


J1/ 3 


-Vx  -W)  F; 


Hi 


(91) 


(92) 


we  rewrite  Eq.  (90)  as 


SuDEu (0)  •  W  =  PirWu. 


(93) 


We  introduce  the  modihed  Cauchy  stress  tensor,  and  push  forward  to  the 
current  conhguration  Eq.  (93) 


(7  = 


-fsft, 


(94) 


ij 
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where 


(95) 


^i,j  j 


1  (i r(J1/3Va'  ■  w) 


Vc  •  w  )  Sij 


(96) 


3  V  J1/ 3 

Note  that,  for  all  geometrical  transformations,  we  use  J  =  detF,  not  detF. 
The  symmetrical  part  of  Eq.  (96)  can  be  written  in  direct  notation: 

(v^w)sym  =  (v*w)‘y”  +  i  _  v-.wf  i, 

=  e(w)  +  ^  O  [e(w)]  -  tr  [e(w)])  I, 

=  e(w)  +  £dll(w)  -  edll(w), 

=  e(w).  (97) 


Thus,  we  can  write  the  stationarity  of  the  modified  potential  energy  given  in 
Eq.  (86)  pushed  forward  to  the  deformed  configuration  as 


0  =  /  (Tij^Eij^dv  -  /  Wifidv  -  /  Wihid'y , 

JB'  JB'  JdhB' 


(98) 


where  /  o  0  =  JF,  and  ( h  o  0)dy  =  TidT.  This  is  the  hnite  deformation  coun¬ 
terpart  of  Eq.  (43),  and  leads  to  the  corresponding  variational  formulation. 


Find  (ft  G  S,  such  that  Vw  e 
a(w,  </>)  =  (w,  f)  +  (w,  h )dhB,, 


(99) 


where  a(w,  <f>)  is  the  first  integral  in  Eq.  (98)  and  V</,  =  {w  o  <p  :  fl 
M3  |  wo  <fi\ dt/>B  =  0}  is  the  tangent  space  to  S  at  <fi . 


4-3  Linearized  operator 


To  solve  the  resultant  nonlinear  problem,  we  derive  the  corresponding  con¬ 
sistent  linearized  operator  (see  Hughes  and  Pister  [36]  and  Sirno  and  Hughes 
[53])  that  will  be  used  in  the  iterative  Newton  algorithm.  Taking  inspiration 
from  the  linearized  operator  for  the  standard  compressible  hnite  deformation 
case  and  the  geometrically  linear  B  theory,  we  might  expect  to  attain  the  usual 
linearized  operator  in  terms  of  “barred”  quantities,  namely, 

/  (WJ  cijid  UkJ  +  Wj  Wj  Ukd)  dv  (100) 

JB' 

However,  as  observed  by  Simo  et  al.  [60] ,  we  also  find  that  the  structure  of  the 
tangent  operator  associated  with  the  projection  method  in  the  fully  nonlinear 
case  leads  to  additional  “unexpected”  terms. 
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To  obtain  the  tangent  operator  we  need  to  compute  the  directional  deriva¬ 
tive  of  the  internal  virtual  work  given  by  the  first  integral  in  Eq.  (98),  or  its 
reference  configuration  version  given  in  Eq.  (90): 


0  [, Su(4 >  +  eU )DEu(<f>  +  eU).W; 


de 


e=0  dV 


(101) 


This  calculation  is  rather  complex  and  is  presented  in  detail  in  Appendix  B. 
The  complete  linearized  operator  is  given  in  the  following  expression: 


IB' 


+  J  <7%j  Q(V*  •  U  -  Vx  •  u )wid  +  ^(V2 


w- V*-  wrn. 


-(V;c  u-Vx  ■  u)(V*  •  w  -  Vx  •  w)Sij 
9 


-(V,T  •  u  Vx  •  w  —  V*  ■  u  •  w  )Sij 


-tr  [V'TuVa’w  -  Va:uVxwJ  Sij  )  dv, 
o 


(102) 


where  the  notations  introduced  in  Eqs.  (B.5),  (B.20)  and  (B.21)  have  been 
used.  Note  that,  contrary  to  the  F  approach  presented  in  de  Souza  Neto  et  al. 
[20,  22],  the  proposed  method  leads  to  a  symmetric  linearized  operator. 


Remark  1  From  the  implementational  point  of  view,  the  same  problems  con¬ 
cerning  the  sparsity  structure  of  the  global  system  that  were  noted  in  the  linear 
case  arise  here  too.  The  “expected”  contributions  to  the  linearized  operator, 
that  is,  the  terms  in  Eq.  (100),  require  a  modification  of  the  sparsity  structure 
of  the  system  equivale?it  to  what  is  needed  in  the  linear  case.  The  additional 
“unexpected”  terms  require  much  more  computational  and  storage  effort  be¬ 
cause  of  terms  involving  projections  of  quantities  that  depend  both  on  the 
weighting  and  solution  spaces.  However,  the  same  strategy  as  in  the  linear  case 
can  be  used.  Although  this  requires  the  assembly  of  more  “mixed”  terms,  the 
computational  and  storage  cost  of  matrix-vector  products  is  in  our  experience 
much  smaller  than  the  computation  and  assembly  of  the  global  matrix. 


The  discrete  form  of  the  consistent  tangent  operator  given  in  Eq.  (102)  is 
presented  in  Appendix  C. 
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5  Numerical  examples  in  the  nonlinear  regime 


The  performance  of  the  proposed  F  methodology  is  numerically  assessed  in  this 
section.  The  first  example,  which  is  a  nonlinear  generalization  of  the  Cook’s 
membrane  problem  presented  earlier,  is  widely  used  in  the  literature.  Other 
two-  and  three-dimensional  examples  found  in  the  literature  are  presented  and 
results  are  compared  with  the  other  techniques  mentioned  in  the  introduction. 


5.1  Cook’s  membrane 


The  setup  of  the  nonlinear  Cooke’s  membrane  example  is  similar  to  the  one 
previously  presented  in  the  geometrically  linear  case.  The  geometry,  loading 
and  boundary  conditions  are  given  in  Figure  8.  This  example  has  been  solved 
in  References  [2,  14,  19,  20,  22,  41,  42,  46,  50,  57]  with  enhanced  assumed 
strain  methods,  the  F  method  of  de  Souza  Neto  et  al,  and  stabilized  mixed 
formulations. 

The  material  model  used  is  Neo-Hookean  following  the  additive  decomposition 
of  the  stored  energy  function  given  in  Eq.  (75).  The  isochoric  and  volumetric 
parts  of  T  are  (see,  e.g.,  [58]): 


T(J,  C)  =  U(J)  +  Tiso(J-2/3C) 

(103) 

u{j)= sKG(j2_i)_inj) 

(104) 

Tiso(J-2/3C)  =  ^/j  ( J~2/3tr[C]  -  3) 

(105) 

where  k  is  the  bulk  modulus  and  /a  the  shear  modulus.  Following  Reference 
[20],  the  values  for  the  material  parameters  are  k  =  40.0942  x  104MPa  and 
/i  =  80.1938M Pa  which  corresponds  to  a  Poisson  ratio  of  u  —  0.4998.  The  load 
value  is  identical  to  the  one  used  in  the  linear  case,  that  is,  F  =  100 N/mm.  The 
quantity  of  interest  is  again  the  vertical  displacement  of  the  top  right  corner 
of  the  plate.  The  convergence  as  a  function  of  the  number  of  elements  per  edge 
and  total  number  of  degrees  of  freedom  are  used  as  criteria  to  investigate  the 
performance  of  the  method.  The  results  for  the  vertical  displacement  versus 
the  number  of  elements  per  edge  are  given  in  Figure  14  and  the  ones  for  the 
vertical  displacement  versus  the  total  number  of  degrees  of  freedom  are  given 
in  Figure  15. 

The  results  are  similar  to  the  ones  obtained  with  the  B  method  in  the  linear 
case.  The  converged  result  obtained  for  quartic  NURBS  functions  is  very  close 
to  values  obtained  with  previously  published  results.  Without  any  treatment  of 
the  incompressibility  condition,  increasing  the  polynomial  order  of  the  NURBS 
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Fig.  14.  Vertical  displacement  of  top  right  corner  versus  number  of  elements  per 
edge  with  and  without  F  for  various  NURBS  orders  obtained  from  ^-refinement. 


Fig.  15.  Vertical  displacement  of  top  right  corner  versus  number  of  degrees  of  free¬ 
dom  with  and  without  F  for  various  NURBS  orders  obtained  from  fc-refinement. 
Note  that  Q\  results  without  F  are  not  shown  on  the  figure. 


37 


functions  improves  the  results,  however,  we  can  see  that  C1  quadratic  NURBS 
still  suffer  from  locking  for  coarse  meshes:  we  need  a  32  x  32  elements  mesh  to 
attain  the  converged  result.  The  use  of  F  significantly  improves  the  results,  and 
good  convergence  is  attained  even  with  Qi/Qo-  Note  that  the  results  obtained 
with  Qi/Qo  are  better  than  Q2  and  that  Q2/Q1  give  better  results  than  Q 3. 
For  high-order  NURBS  with  F,  we  obtain  converged  results  with  very  coarse 
meshes.  The  results  in  term  of  displacement  versus  total  number  of  degrees 
of  freedom  show  again  the  advantages  of  ^-refinement.  Note  that  in  terms  of 
unknowns,  Q2/Q1  is  more  accurate  than  Q4.  Qi/Q2  and  Qa/Qz  give  better 
results  still,  but  the  improvement  over  Q2/Q\  is  not  as  great  as  in  the  linear 
case. 


Iteration  number 

Relative  norm  of  the  residual 

1 

l.OOOOOOOOxlO-1 

2 

9.99179792x  10-2 

3 

1.31723792x  10-2 

4 

9.23073640x  10~3 

5 

3.33352491  xl0“7 

Table  2 

Cook’s  membrane.  Evolution  of  the  relative  norm  of  the  residual  during  the  last 
Newton  step  for  a  8  x  8  elements  mesh  with  the  proposed  F  technique  with  Qi/Qo- 

We  finally  investigate  the  behavior  of  the  consistently  linearized  operator  de¬ 
scribed  previously.  We  see  in  Table  2  the  evolution  of  the  relative  Euclidian 
norm  of  the  residual  over  the  iterations  for  the  last  loading  step  for  an  8x8 
element  mesh  with  the  proposed  F  technique  with  Q\/Qq.  This  illustrates  that 
the  tangent  matrix  obtained  from  the  consistent  linearization  of  the  problem 
results  in  quadratic  convergence  rate.  The  doubling  of  zeros  behind  the  decimal 
point  as  we  get  very  close  to  the  solution  is  typically  seen  in  our  calculations. 


5.2  Plane  strain  compression  of  a  block 


The  next  example  is  the  plane  strain  compression  of  a  rectangular  block.  This 
problem  has  been  studied  in  Reference  [50]  as  a  severe  benchmark  for  enhanced 
element  formulations  as  they  often  exhibit  locking  in  such  situations  [2,  21], 
The  geometry,  loading  and  boundary  conditions  are  given  in  Figure  16.  For 
symmetry  reasons,  only  one  half  of  the  specimen  needs  to  be  considered. 

The  material  model  is  Neo-Hookean  with  a  slightly  different  stored  energy 
function  than  the  one  used  in  the  previous  example: 

T(  J,  C)  =  ^/i  (tr[C]  —  3)  —  n  In  J  +  ( J2  -  1  -  2  In  J )  (106) 
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Fig.  16.  2D  plane  strain  compression  test:  geometry,  boundary  conditions,  loading, 
and  typical  results  of  Von  Mises  Stress  field  on  the  deformed  configuration. 


Following  Reference  [50],  the  values  for  the  material  parameters  are  k  = 
400889.806 MPa  and  /i  =  80.194MPa.  The  load  factor  p/p0  is  increased  by 
steps  of  0.5  with  a  reference  value  of  po  =  20 MPa  (a  unit  thickness  is  consid¬ 
ered).  Note  that  a  “dead”  load  is  considered  as  shown  in  Figure  16,  that  is, 
one  assumes  a  fixed  vertical  surface  load  specified  on  the  reference  configura¬ 
tion  during  the  entire  deformation  of  the  structure.  The  Dirichlct  boundary 
conditions  are  the  following:  no  vertical  displacement  on  the  bottom  and  no 
horizontal  displacement  on  the  top.  The  quantity  of  interest  is  the  compres¬ 
sion  level  of  the  upper  middle  point.  The  evolution  of  the  compression  level 
versus  the  number  of  elements  per  edge  for  loading  ratios  from  p/po  =  10  to 
p/Po  —  60  is  studied. 


In  Figure  17,  we  compare  the  results  obtained  with  the  proposed  F  approach 
with  linear  elements  (he.,  Qi/Qo)  with  the  mixed  element,  referred  to  as  Qi/Po 
(on  the  left),  and  with  the  enhanced  element  QiSP,  both  from  [50].  For  all 
loading  levels,  convergence  is  attained  for  meshes  of  16  x  16  elements  or  less. 
Note  that,  with  coarse  meshes,  the  F  method  performs  better  than  both  the 
Qi/Po  and  QiSP  elements  at  high  compression  levels  while  at  low  compression 
levels  the  situation  is  reversed.  Results  obtained  with  the  proposed  F  approach 
for  linear,  quadratic  and  cubic  NURBS  (i.e.,  Qi/Qo ,  Q2/Q1  and  Q3/Q2)  are 
then  compared  in  Figure  18.  Although  quadratic  and  cubic  functions  tend  to 
produce  a  stiffer  response  than  linears,  we  can  see  that  in  all  cases  we  reach 
convergence  for  16  x  16  element  meshes.  Quadratic  and  cubic  NURBS  produce 
better  results  for  coarser  meshes.  Note  that  the  jump  from  Q2/Q1  to  Q3/Q2 
only  slightly  improves  the  results  as  was  observed  for  the  previous  example. 
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Fig.  17.  Compression  level  in  %  versus  number  of  elements  per  side  for  the  two-di¬ 
mensional  plane  strain  block  for  various  load  levels.  Left  side:  F  Qi/Qo  compared  to 
the  mixed  u-p  formulation  Qi/Pq  (results  from  [50]).  Right  side:  F  Q\/Qq  compared 
to  the  stabilized  element  Q\SP  (results  from  [50]). 

5.3  Three-dimensional  compression  of  a  block 


The  next  example  is  the  three-dimensional  generalization  of  the  previous  one. 
It  was  studied  in  Reference  [51]  as  a  benchmark  for  the  three-dimensional  ver¬ 
sion  of  the  QiSP  element.  This  problem  is  of  interest  to  test  the  performance 
of  enhanced  elements  in  three  dimensional  compression  cases.  The  geometry, 
loading  and  boundary  conditions  are  given  in  Figure  19,  along  with  a  typi¬ 
cal  plot  of  the  Von  Mises  stress  field  on  the  final  deformed  configuration.  For 
symmetry  reasons,  only  one  quarter  of  the  specimen  is  considered. 

The  material  model  is  Neo-Hookean  with  a  somewhat  different  stored  energy 
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Fig.  18.  Compression  level  in  %  versus  number  of  elements  per  side  for  the  two-di¬ 
mensional  plane  strain  block  for  various  load  levels.  Comparison  of  F  results  for 
various  NURBS  orders  obtained  from  ^-refinement. 


Fig.  19.  3D  compression  test:  geometry,  boundary  conditions,  loading,  and  typical 
results  of  the  Von  Mises  stress  field  on  the  deformed  configuration. 


function  than  the  previous  ones: 


*(4C) 


^/i  (tr[C]  —  3)  —  /i  In  J  +  (In  J)“ 


(107) 


Following  Reference  [51]  the  values  for  the  material  parameters  are  k  = 
400889. 806MPa  and  /i  =  80.194MPa.  Loading  factors  from  p/p0  =  20  to 
p/po  =  80  are  studied  with  a  reference  value  of  p0  =  AMPa.  Note  that  a 
“dead”  load  is  considered  as  shown  in  Figure  19,  that  is,  one  assumes  a  fixed 
vertical  surface  load  specified  on  the  reference  configuration  during  the  entire 
deformation  of  the  structure.  The  Dirichlet  boundary  conditions  are:  no  ver- 
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tical  displacement  on  the  bottom  surface  and  no  horizontal  displacement  on 
the  top  surface. 
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Fig.  20.  Compression  level  in  %  versus  number  of  elements  per  side  for  the  three-di¬ 
mensional  block  for  various  load  levels.  Left  side:  F  Qi/Qq  compared  to  the  mixed 
u-p  formulation  Qi/Pq  ([51]).  Right  side:  F  Qi/Qo  compared  to  the  stabilized  ele¬ 
ment  QiSP  ([51]). 

The  quantity  of  interest  is  the  compression  level  of  the  upper  center  point.  The 
evolution  of  the  compression  level  versus  the  number  of  elements  per  edge  for 
various  loading  ratios  is  studied. 

In  Figure  20,  we  compare  the  results  obtained  with  the  proposed  F  approach 
with  linear  elements  (he.,  Qi/Qo)  with  the  mixed  element  Q \ / Pq  (on  the  left) 
and  with  the  enhanced  element  QiSP,  both  from  [51].  We  can  see  that  con¬ 
vergence  is  obtained  for  meshes  of  8  x  8  x  8  elements  and  that  with  coarser 
meshes  the  results  obtained  with  Qi/Qo  are  almost  identical  to  the  ones  ob¬ 
tained  with  Q\/Pq.  Results  obtained  with  linear,  quadratic  and  cubic  NURBS 
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Fig.  21.  Compression  level  in  %  versus  number  of  elements  per  side  for  the  three-di¬ 
mensional  block  for  various  load  levels.  Comparison  of  F  results  for  various  NURBS 
orders  obtained  from  /c-refinement. 

are  presented  in  Figure  21.  As  for  the  two-dimensional  case,  quadratic  and 
cubic  functions  tend  to  give  a  stiffer  response  than  linear  functions.  Never¬ 
theless,  results  are  sligthly  better  with  the  higher-order  functions  than  with 
linear  functions.  The  improvement  obtained  by  going  from  quadratic  to  cubic 
NURBS  is  quite  pronounced,  especially  for  the  coarse  2x2x2  element  mesh. 


5-4  Thick  cylindrical  shell 


I11  the  fourth  example,  we  consider  a  geometry  which  cannot  be  represented 
exactly  by  piecewise  linear  polynomials.  It  involves  the  compression  of  a  thick 
hollow  cylinder,  for  which  the  exact  initial  geometry  can  be  obtained  with 
quadratic  NURBS.  This  example  was  treated  in  [13,  51]  and  interest  was 
focused  on  comparisons  with  shell  formulations  and  shear  locking.  As  we 
have  seen  in  the  Cooke’s  membrane  problem,  the  proposed  method  is  able 
to  avoid  both  volumetric  and  shear  locking.  We  note  for  this  case  that  the 
usual  Galerkin  formulation  using  cubic  and  quartic  functions  also  produced 
good  results,  without  evidence  of  volumetric  locking.  The  reason  for  this  is 
that  the  traction-free  cylindrical  surfaces  alleviate  the  tendency  to  lock.  This 
tendency  is  also  seen  in  plane  stress  analyses  where  there  is  no  locking  problem 
whatsoever.  However,  using  the  F  method  in  conjunction  with  these  high-order 
functions  improves  the  results  further.  Therefore,  we  consider  the  thick  cylin¬ 
der  with  a  compressible  material  in  a  bending  case  where  shear  locking  is 
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Fig.  22.  Thick  cylindrical  shell:  geometry,  loading,  boundary  conditions  and  typical 
results  for  the  shear  component  ayz  of  the  Cauchy  stress  tensor  on  the  deformed 
configuration  (note  that  the  scale  has  been  reduced  to  see  the  interesting  features). 

expected  to  occur  to  further  assess  this  observation. 

The  geometry,  boundary  conditions,  loading,  and  a  typical  result  of  the  shear 
stress  component  ayz  of  the  Cauchy  stress  tensor  on  the  final  deformed  config¬ 
uration  are  given  in  Figure  22.  As  in  [51],  we  use  only  one  quadratic  element 
through  the  thickness,  which  is  enough  to  capture  the  bending  behavior.  Be¬ 
cause  of  symmetry,  only  one  fourth  of  the  structure  is  considered.  The  mate¬ 
rial  model  considered  is  the  same  as  in  the  previous  example  with  parameters 
k  =  240000MPa  and  //  =  6000MPa,  which  corresponds  to  a  Poisson’s  ratio 
of  v  =  0.4,  therefore  no  volumetric  locking  is  expected.  Note  that  again  a  dead 
load  is  considered  as  shown  in  Figure  22,  that  is,  one  considers  a  fixed  vertical 
line  load  specified  on  the  reference  configuration  during  the  entire  deformation 
of  the  structure.  The  quantity  of  interest  is  the  vertical  displacement  of  point 
A  shown  in  Figure  22. 

The  results  of  the  vertical  displacement  at  point  A  are  shown  in  Figure  23, 
where  we  compare  the  results  obtained  with  and  without  F  with  C1-quadratic 
and  C2-cubic  NURBS,  Q\SP  and  the  enhanced  assumed  strain  nonlinear  shell 
element  proposed  in  [13],  versus  the  number  of  elements  in  the  circumferential 
direction.  As  in  [51],  we  used  the  following  sequence  of  meshes  (circumference 
x  thickness  x  axis):  4 x  1  x  2,  8  x  1  x  4,  16  x  1  x  8,  32  x  1  x  16  that  were  obtained 
from  ^-refinement.  We  can  see  that  all  cases  are  converging  to  approximately 
the  same  value  obtained  in  [51]  with  both  Q\SP  and  shell  elements,  and  that 
cubic  functions  give  better  results  especially  for  the  coarsest  meshes.  The  use 
of  F  with  quadratic  functions  improves  the  results  by  approximatively  50%  for 
the  coarse  mesh,  and  10%  for  the  next  mesh.  Note  also  that,  for  cubic  NURBS, 
the  use  of  F  improves  the  results  for  the  coarse  mesh  by  approximatively  10%, 
and  that  the  two  curves  collapse  for  the  third  discretization. 
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Fig.  23.  Vertical  displacement  of  point  A  versus  number  of  elements  in  the  circum¬ 
ferential  direction  with  and  without  F  for  various  NURBS  orders  obtained  from 
^-refinement,  for  Q±SP  (results  from  [51])  and  for  an  enhanced  assumed  strain 
nonlinear  shell  element  (results  from  [13,  51]). 

5.5  Pinched  torus 


The  last  example  again  exploits  the  ability  of  NURBS  to  exactly  represent 
conic  sections.  It  consists  of  the  pinching  of  a  toroidal  solid,  and  is  similar  to 
an  example  proposed  in  [14]. 


Shear  modulus  p 

5.67  MPa 

Bulk  modulus  k 

2.8333  103  MPa 

Inner  radius  r 

8  m 

Outer  radius  R 

10  m 

Reference  pressure  po 

0.195  MPa 

BC  in  plane  x  =  0 

ux  =  0 

BC  in  plane  y  =  0 

o 

II 

53 

BC  in  plane  z  =  0 

O 

II 

3 

Table  3 

Pinched  torus:  material  properties  and  boundary  conditions. 

The  geometry,  loading,  boundary  conditions  and  mesh  are  shown  in  Figure  24, 
and  the  material  parameters  are  given  in  Table  3.  The  same  Neo-Hookean  ma¬ 
terial  model  as  for  the  Cook’s  membrane  is  considered,  which  corresponds  to  a 
Poisson’s  ratio  of  u  —  0.4998.  Due  to  symmetry  condition,  only  one  quarter  of 
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Fig.  24.  Pinched  torus:  geometry,  quarter  mesh,  loading  and  boundary  conditions. 

the  structure  is  considered,  with  the  corresponding  symmetry  boundary  con¬ 
ditions  applied.  The  quarter  mesh  with  4  x  16  x  2  elements  (that  is  4  elements 
in  the  “large”  circumferential  direction,  16  elements  in  the  “small”  circum¬ 
ferential  direction  and  2  elements  in  the  radial  direction)  shown  in  Figure  24 
with  quadratic  and  cubic  NURBS  is  used. 

The  axx  and  azz  components  of  the  Cauchy  stress  tensor  plotted  on  the  fi¬ 
nal  deformed  configuration  with  and  without  F  for  C1-quadratic  and  C2-cubic 
functions  are  shown  in  Figure  25.  The  differences  on  the  final  deformed  con¬ 
figurations  show  that  Qi  without  F  suffers  form  locking.  Although  F  Q2/Q1 
and  Q3  without  F  look  similar,  we  can  see  that  F  Q3/Q2  improves  the  result. 
The  stress  contours  for  the  component  crxx  of  the  Cauchy  stress  tensor  shows 
typical  oscillations  due  to  locking  for  Qo  without  F.  These  oscillations  are  not 
observed  in  the  three  other  cases  for  this  component.  However,  we  can  see  on 
the  bottom  part  of  Figure  25,  that  oscillations  are  present  for  both  quadratic 
and  cubic  meshes  without  F  for  the  component  azz  of  the  Cauchy  stress  ten¬ 
sor.  Note  that,  for  both  components  and  both  orders  of  approximation,  the 
results  with  F  do  not  present  such  oscillations:  the  stresses  are  smooth  and 
very  similar  for  Q2/Q1  and  Q‘$IQ2- 

The  vertical  displacement  at  point  A  for  quadratic  and  cubic  functions  with 
and  without  F  versus  the  loading  step  is  shown  in  Figure  26.  We  also  see  that 
Q2  suffers  from  locking  and  that  at  the  final  load  the  difference  between  Q2  and 
Q2/Q1  is  approximative^  of  50%.  As  we  increase  the  order  of  approximation 
without  F,  we  see  that  the  results  converge  to  those  of  Q2/Q i-  Using  C 2  cubic 
functions  without  F  significantly  improve  the  results  which  are  very  close  to 
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Fig.  25.  Pinched  torus:  component  (a)  axx  and  (b)  azz  of  the  Cauchy  stress  tensor 
on  the  deformed  configuration  with  and  without  F  for  (^-quadratic  functions  and 
C2- cubic  functions  (note  that  the  scales  have  been  reduced  to  see  the  interesting 
features). 
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those  of  Q2/Qi  -  However,  we  can  see  that  F  Q3IQ2  improves  the  result  further. 
Note  that,  /i-refining  the  mesh  in  the  circumferential  direction,  that  is,  using 
a  mesh  of  8  x  16  x  2  elements,  and  F  Q2/Q1  allows  us  to  obtain  exactly  the 
same  results  as  for  F  Q3/Q2  on  the  cruder  mesh. 


Fig.  26.  Pinched  torus:  vertical  displacement  at  point  A  versus  load  step  number 
with  and  without  F  for  various  NURBS  orders  obtained  from  fc-refinement. 


5.6  Elastic-plastic  Cook’s  membrane 


We  present  some  preliminary  results  for  elastic-plastic  behavior  at  finite  strain 
with  the  proposed  F  approach.  We  consider  the  same  hyperelastic  extension  of 
J2-flow  theory  with  poly-convex  hyperelastic  uncoupled  stored  energy  function 
and  nonlinear  isotropic  hardening  as  in  Simo  [55,  56]  and  Simo  and  Hughes 
[58].  It  consists  of  a  Neo-Hookean  model  for  the  elastic  part,  and  an  associative 
flow  rule  based  on  a  Von  Mises  yield  criterion  with  isotropic  hardening  follow¬ 
ing  a  saturation  law  for  the  plastic  part.  The  nonlinear  isotropic  hardening 
rule  is  defined  by  the  following  equation: 

k(a)  =  a0  +  (doc  —  cr0)[l  —  exp(— 5a)]  +  Ka,  with  5  >  0,  (108) 

where  <r0  is  the  initial  flow  stress,  is  the  saturation  flow  stress,  5  is  the  sat¬ 
uration  exponent,  K  is  the  linear  hardening  coefficient,  and  a  is  the  equivalent 
plastic  strain. 
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The  Neo-Hookean  material  model  and  boundary  conditions  are  identical  to 
the  ones  previously  used  for  this  problem,  and  the  load  value  is  chosen  to  be 
F  =  5N/mm.  The  material  parameters  are  those  used  in  Simo  and  Armero 
[57],  and  are  presented  in  Table  4. 


Shear  modulus  /jl 
Bulk  modulus  k 
Initial  flow  stress  <7o 
Saturation  flow  stress  aQ 0 
Saturation  exponent  5 
Linear  hardening  coefficient  I\ 


80.1938  MPa 
164.21  MPa 
0.450  MPa 
0.715  MPa 
16.93 

0.12924  MPa 


Table  4 

Elastic-plastic  Cook’s  membrane:  material  properties. 


The  results,  shown  in  Figure  27  for  the  vertical  displacement  of  the  top  right 
point  versus  number  of  elements  per  side,  compare  well  to  those  obtained  in  [2, 
57] .  Similar  convergence  behavior,  when  increasing  the  order  of  approximation 
with  the  proposed  F  method,  to  those  obtained  in  the  nearly  incompressible 
linear  and  nonlinear  clastic  cases  is  obtained.  One  can  also  see  the  equivalent 
plastic  strain  for  a  16  x  16  element  mesh  with  F  using  cubic  NURBS  functions 
in  Figure  28. 
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Fig.  27.  Nonlinear  elastic-plastic  Cook’s  membrane  with  nonlinear  isotropic  hard¬ 
ening.  Vertical  displacement  of  top  right  corner  versus  number  of  elements  per  edge 
with  and  without  F  for  various  NURBS  orders  obtained  from  fc-refinement. 
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Fig.  28.  Nonlinear  elastic-plastic  Cook’s  membrane  with  nonlinear  isotropic  hard¬ 
ening.  Final  equivalent  plastic  strain  with  F  for  a  mesh  of  16  x  16  elements  using 

Q‘i/Q2- 


6  Conclusions 


In  this  paper  we  have  presented  B  and  F  projection  methods  for  nearly  incom¬ 
pressible  small  deformation  and  large  deformation  elasticity  in  the  framework 
of  isogeometric  analysis.  These  methods  make  use  of  the  higher  continuity 
property  of  high-order  NURBS  obtained  from  ^-refinement,  a  unique  feature 
of  isogeometric  analysis.  We  first  presented  an  extension  of  the  well  known 
B  method  in  this  context.  This  allowed  us  to  concentrate  on  the  choice  of 
projection  spaces  and  techniques  to  solve  problems  with  higher-order  func¬ 
tions.  The  results  obtained  showed  the  efficacity  of  the  projection  methods 
to  solve  nearly  incompressible  problems  with  a  pure  displacement  formulation 
using  higher-order  functions  with  higher  continuity.  The  small  increase  in  the 
number  of  unknowns  when  increasing  the  order  of  approximation  within  the 
fc- refinement  process  (contrary  to  p- refined  finite  elements)  was  also  shown  to 
be  an  interesting  feature  of  isogeometric  analysis. 

This  study  demonstrates  that  the  good  properties  of  NURBS  previously  ob¬ 
served  in  other  applications  are  also  present  with  projection  methods.  The 
recent  work  on  the  use  of  high  continuity  functions  (especially  C 1  quadrat¬ 
ics)  compared  with  C°  functions,  in  fluid  mechanics  [1,  4,  5]  and  structural 
vibrations  [18],  gives  us  a  better  understanding  of  the  effect  of  increased  con¬ 
tinuity  on  approximation.  We  anticipate  that  such  properties  could  provide 
interesting  improvements  in  large  scale  nonlinear  solid  mechanics. 
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As  a  first  step,  we  proposed  a  new  F  projection  method  to  treat  the  in¬ 
compressibility  constraint  in  the  geometrically  nonlinear  regime.  Although  it 
shares  features  with  the  F  method  of  de  Souza  Neto  et  al.  [20,  22],  we  think 
of  it  as  a  more  theoretically  sound  method  since  it  is  based  on  a  modified 
minimum  potential  energy  principle  and  leads  to  a  symmetric  consistent  tan¬ 
gent  which  is  not  the  case  with  the  method  of  de  Souza  Neto  et  al.  Moreover, 
the  method  does  not  make  any  assumption  on  the  polynomial  order  of  the 
approximation  space  and  is  therefore  not  limited  to  piecewise  linear  elements. 

The  numerical  results  show  that  the  method  compares  well  with  previously 
published  methods,  especially  for  the  Cook’s  membrane  problem  for  which  we 
obtained  virtually  identical  results  to  those  of  [20,  22]  with  piecewise  linear 
functions.  We  also  investigated  two-  and  three-dimensional  compression  prob¬ 
lems  where  the  limits  of  enhanced  strain  elements  are  usually  reached,  and 
we  observed  good  performance  of  the  F  isogeometric  formulation  even  with 
low  order  functions.  The  ability  of  the  method  to  overcome  shear  locking  for 
quadratic  and  higher-order  basis  functions  was  also  demonstrated  although 
no  particular  attention  was  paid  to  that  topic  in  the  design  of  the  method. 
Finally,  some  preliminary  calculations  considering  an  elastic-plastic  behavior 
coupled  with  the  proposed  F  method  were  presented,  and  show  very  promising 
results. 

The  proposed  approach  takes  advantage  of  a  simple  pure  displacement,  pro¬ 
jection  method  and  high-order  fc-rehned  NURBS  in  the  isogeometric  analysis 
framework.  It  opens  a  new  door  to  the  use  of  precise  CAD  geometry  rep¬ 
resentations  and  higher-continuity  properties  of  NURBS-based  isogeometric 
analysis  in  nonlinear  structural  mechanics  applications.  Areas  of  further  re¬ 
search  in  structural  mechanics  needed  include  continuing  the  investigation  of 
the  approach  initiated  herein  with  clastic-plastic  materials,  and  extending  it 
to  frictional  contact. 
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A  Directional  derivative  of  the  modified  Green-Lagrange  strain 
tensor 


To  obtain  the  directional  derivative  of  E,  we  first  need  to  compute  the  deriva¬ 
tive  of  the  modified  deformation  gradient  F : 


dF (0  +  eW)  _  d 
Fe  |e=0  "  Fe 


JV  3(0 +  eW) 
JV  3(0 +  eW) 


F(0  +  eW) 


e=0 


(A.l) 


Using  the  linearity  of  the  projection  operator  and  expanding  the  former  ex¬ 
pression,  we  get: 


d¥{<j>  +  eW) 

Fe 


e=0 


3JV3(0)  1 

JV3(0) 
+  J1/ '3(0) 


7T  •  W)  -  J1/3(0)Va; 

VAW 


•  w 


F 

(A.2) 


For  simplicity,  we  make  use  of  the  notation  introduced  in  Eq.  (73): 


J1/  3(0) 
J1/3(0)1 


(A.3) 


and  therefore  Eq.  (A.2)  becomes: 


DF(0)-W  =  aVAW  +  - 

3 


,  1  ( vr(J1/3(0)V-  ■  W)  „ 


J1/ 3(0) 


—  Vx  •  W  F. 


(A.4) 


The  definition  of  the  modified  Green-Lagrange  strain  tensor  in  terms  of  the 
modified  deformation  gradient  is  used  in  conjunction  with  the  previous  equa¬ 
tion: 


2DE(0)  •  W  =  (DF(<f>)  ■  W)T F(0)  +  FT(0)DF(0)  •  W,  (A.5) 
which  finally  allows  us  to  write: 

DE(<f>)  •  W  =  -  ( 7r^J  l—*  _  yT  ■  W^)  FtF  +  a  (FTVA'W)sym 

3  \  J V3(0)  / 

(A.6) 
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B  Consistent  linearized  operator 


We  begin  by  recalling  Eq.  (101): 


d  [Su{<f>  +  eU )DEu(tf>  +  eU)  •  W 


.)L.U  {DEkl  •  U)  (DEjj  •  w) 

+  SjjD  (FtJW-r)  (B.l) 


We  focus  on  the  first  term  of  this  equation.  We  define  the  modified  fourth 
order  material  tensor  as  we  previously  defined  Eq.  (89),  the  modified  second 
Piola-Kirchhoff  stress  tensor: 
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(B.2) 


Using  the  result  from  Appendix  A  (he.,  Eq.  (A. 6)),  the  first  term  of  Eq.  (B.l) 
becomes 

C  ijkl  F%j  Wtj  Fj  tUjK.  (B.3) 

As  done  previously,  we  can  push  forward  this  equation  to  the  current  config¬ 
uration  using: 


Wij  =  wi)kFkI  =  a  (wi,k  +  ^  [V* 


w 


w]  Sik 


kl  i 


(B.4) 


where  we  have  introduced  for  brevity  the  following  notation  in  both  the  ref¬ 
erence  and  current  configuration: 


— —  7r(./1/3(0)V'r  •  W)  7T (JVS^V*  •  w)  - - 

V*  •  W  =  — — === - -  =  — — ==A - -  =  Vx  ■  w  (B.5) 

J1/ 3(0)  JV3(0) 


The  push  forward  of  Eq.  (B.3)  is: 

C ' i jklFu FjjFkK FuM j  UkJ.  (B.6) 


We  can  consider  the  classical  push  forward  on  the  current  configuration  of  the 
modified  material  tensor  in  terms  of  “bar”  quantities,  that  is: 


Cijki  —  ~j  C  j  j  k  l  Ft  i  Fj  j  Fk  kFiLi 


(b.t) 


and  finally  obtain  the  expected  material  contribution  to  the  linearized  operator 
on  the  current  configuration: 


U'  ®  J^ijM  k  jd  v 


(B.8) 
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We  now  consider  the  second  term  in  Eq.  (B.l),  which  can  in  turn  be  split  into 
two  terms: 


SjjD  (FijWij)  •  U  =  Sjj  [DFU  ■  U]  WtJ  +  SuFtJ  [DWh,  ■  U]  (B.9) 

The  first  term  represents  the  geometrical  contribution  to  the  linearized  oper¬ 
ator  in  terms  of  “bar”  quantities,  and  the  second  term  is  the  one  responsible 
for  the  additional  contribution.  We  first  focus  on  the  geometrical  contribution. 
Using  the  result  obtained  in  Appendix  A  (he.,  Eq.  (A. 6)),  we  have: 

DFiJ-XJ  =  Uij,  (B.10) 

which  can  be  introduced  in  Eq.  (B.9),  using  the  definition  of  the  modified 
Cauchy  stress  tensor  given  in  Eq.  (94): 

b'/  jl't.j  !  U.  /  JcTijWkj  Uk,i-  (B.ll) 

The  geometrical  contribution  to  the  linearized  operator  finally  takes  the  ex¬ 
pected  form: 

/  wi~GijU]~dv.  (B.12) 

JB' 

To  determine  the  second  term  in  Eq.  (B.9),  we  need  to  compute: 

DW~i-  U,  (B.13) 

where 

W~J  =  l-(W^W  -  V*  •  W )FUFU  +  aFuWu  (B.14) 

Rewriting  the  divergence  of  W  as 

V*  •  W  =  witi  =  WijF-1  =  tr(VA'W  F"1),  (B.15) 

we  get 

D(VX  •  W)  •  U  =  tr(VAW  •  D{ F^1))  U.  (B.16) 

Since  we  have  the  following  equation  for  the  directional  derivative  of  the  inverse 
of  the  deformation  gradient: 

D{ F^1)  •  U  =  — F_1  VA'U  F^1,  (B.17) 

Eq.  (B.16)  can  be  simplified  as: 

D{VX  •  W)  •  U  =  -tr  [V'rwVau] .  (B.18) 

Proceeding  analogously,  we  can  determine  the  next  expression: 

d(WW)  •  U  =  D  ^ J  /-Z_’  —  ^  •  U.  (B.19) 
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Following  Eq.  (B.5),  we  introduce  the  following  notations: 
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which  allows  us  to  express  Eq.  (B.19)  as: 
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We  can  gather  all  terms  to  obtain  the  following  expression  for  Eq.  (B.13): 


DWtj  ■  U  =  -(Vx  ■  u  -  Vx  •  u)V%  F  +  -(Vx  •  w  -  V*  •  w)V:cu  F 
3  3 


-(Vx  •  u  -  V*  •  u)(Vx  •  w  -  Vx  •  w)F 
9 

1  / 1 


(B.23) 
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3  V  3  3 


Using  this  expression  we  finally  obtain  the  last  unexpected  contribution  to  the 
linearized  operator: 


°ij  Q(VX  •  u  -  V'T  •  u)witj  +  I(V- 
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-tr  [V^uV^w  -  VxuVxwJ  Sij  )  dv 
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C  Discrete  consistent  linearized  operator 


To  express  the  discrete  version  of  the  consistent  linearized  operator  given  in 
Eq.  (102),  we  start  by  introducing  the  following: 

Cujj  =  CKILJFiKFjL.  (C.l) 
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The  “material”  contribution  to  the  discrete  consistent  linearized  operator  (i.e., 
the  first  term  in  Eq.  (102))  is 
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The  “geometrical”  contribution  to  the  discrete  consistent  linearized  operator 
(i.e.,  the  second  term  in  Eq.  (102)),  and  parts  of  the  “unexpected”  contribu¬ 
tions  are: 
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Finally,  the  remaining  terms  of  the  “unexpected”  contribution  are: 
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The  corresponding  discrete  internal  forces  can  be  expressed  as: 
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